<figure><img src="image_1.png" /></figure>

HAL
open science

# Robust estimation of regional climate change: construction of an hybrid approach between deep neural networks and climate models

Antoine Doury

► To cite this version:

Antoine Doury. Robust estimation of regional climate change: construction of an hybrid approach between deep neural networks and climate models. Climatology. Université Toulouse 1 Capitole (UT1 Capitole), 2022. English. NNT: tel-04239694

HAL Id: tel-04239694

https://hal.science/tel-04239694

Submitted on 12 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_2.png" /></figure>

Université
de Toulouse

THÈSE

En vue de l'obtention du

DOCTORAT DE L'UNIVERSITÉ DE
TOULOUSE

Délivré par : l'Université Toulouse 1 Capitole (UT1 Capitole)

Présentée et soutenue le 16/12/2022 par :

Antoine DOURY

Robust estimation of regional climate change: construction
of an hybrid approach between deep neural networks and
climate models.

JURY

MATHIEU VRAC

Rapporteur

PATRICK GALLINARI

Rapporteur

ERIKA COPPOLA

Examinatrice

ERWAN LEPENNEC

Examinateur

JOSÉ MANUEL

Examinateur

GUTIÉRREZ

SEBASTIEN GADAT

Directeur de These

SAMUEL SOMOT

Co-Directeur de These

École doctorale et spécialité :

MITT : Domaine Mathématiques : Mathématiques appliquées

Unité de Recherche :

Centre National de Recherche Météorologiques (UMR 3589)

Directeur(s) de Thèse :

Sebastien GADAT et Samuel SOMOT

Rapporteurs :

Mathieu VRAC et Patrick GALLINARI
Abstract

An essential challenge for the climate science community is to provide trustful
information about the local impacts of global warming. Climate models are the
main tool to study climate evolution according to human activities and greenhouse
gas emissions scenarios. They are a numerical representation of the Earth System.
Global climate models (GCM) produce worldwide simulations at too low resolution
to correctly represent extreme meteorological events that strongly impact our so-
cieties. Today we use dedicated regional climate models (RCM) that transform a
global low-resolution simulation into a high-resolution one over an area of interest.
Nevertheless, the high resolution of those models implies a (much) higher cost that
strongly limits the number of those climate simulations and, thus, the necessary
exploration of the different sources of uncertainties.

The present work aims to propose a strategy to recreate, at low-cost, high-resolution simulations from low-resolution ones. The RCM-emulator introduced here aims to estimate the downscaling function included in a RCM using the recent development of neural networks. This study introduces the concept of RCM-emulator and presents a framework to build, train, and evaluate it. The main result of this study is that the RCM emulator is a credible approach to take up this challenge. Indeed it shows an excellent ability to create realistic high-resolution temperature and precipitation fields, consistent with the low-resolution simulation it downscales. We also study the applicability of the tool to the various low-resolution simulations. Moreover, this work also highlights the decisive advantage of using RCM simulations to *learn* this relationship as it allows to explore future climates and poorly known regions.

The conclusions of this study open the door to further development and various
promising applications. Indeed, the RCM-Emulator makes possible the production
of robust messages about the local impacts of climate change. Moreover, another
significant result of this work is that the emulator performance relies strongly on the
calibration set. It is then essential to design the best simulation set to have the most
robust emulator implying maybe to revisit the way of choosing which simulation to
make with a RCM.
# Résumé

Un des grands défis pour la communauté scientifique du climat est de produire une information fiable concernant les impacts locaux du changement climatique global. Les modèles de climat sont les principaux outils pour étudier son évolution en fonction de scénarios liés aux activités humaines et émissions de gaz à effet de serre. Les modèles de climat globaux (GCM) produisent des simulations couvrant l'ensemble du globe à des résolutions trop faibles pour représenter correctement certains évènements météorologiques extrêmes qui impactent fortement nos sociétés. Aujourd'hui, nous utilisons des modèles de climat régionaux (RCM) pour transformer une simulation globale à basse résolution en simulation à haute résolution sur une zone géographique d'intérêt. Cependant, la haute résolution de ces modèles implique des coûts de calcul très élevés limitant le nombre de simulations réalisables et donc l'exploration nécessaire des différentes sources d'incertitudes.

Le présent travail a pour objectif de proposer une stratégie pour créer à moindre coût, des simulations haute résolution à partir d'autres à basse résolution. Le RCM-emulateur, introduit ici, cherche à estimer la fonction de descente d'échelle incluse dans un RCM en utilisant notamment les récents développement concernant les réseaux de neurones. Cette étude introduit le concept de l'émulateur et présente un cadre pour le construire, l'entraîner et l'évaluer. Le résultat principal de cette étude est que l'émulateur est une approche crédible pour relever ce défi. En effet, il montre une excellente capacité à créer des champs de température et précipitation à haute résolution cohérent avec la simulation basse résolution d'origine. Nous étudions aussi son applicabilité à différentes simulation basse résolution. Ce travail met aussi en évidence l'avantage décisif d'utiliser des simulations RCM pour apprendre cette relation puisque cela permet d'explorer des climats futurs et des régions mal observées.

Les résultats de cette étude ouvrent la porte à de nouveau développements mais aussi à diverses applications prometteuses. En effet, le RCM-emulateur rend possible la production de messages robustes concernant les impacts locaux du changement climatique. De plus, un autre résultat important de ce travail est l'importance du jeu de données d'apprentissage sur les performances de l'émulateur. Il est essentiel de définir le meilleur ensemble de simulation pour construire un émulateur robuste ce qui implique peut-être de revoir la manière de choisir quelle simulation faire avec un RCM.
Acknowledgements

I would like to apologize any non-french reader as the aknowledgements will be
in french, it was complicated enough like that... Je n'aurai jamais pensé que cette
page eut été la plus difficile à écrire de cette thèse. En fait si je dois retenir une chose
de ces années, c'est certainement la richesse des rencontres qui ont toutes contribuées
à faire de moi qui je suis et donc de cette thèse ce qu'elle est. Ces remerciements
sont un exercice difficile car à l'opposé du reste de ce manuscrit, alors on va essayer
de faire les choses bien et dans l'ordre.

Je souaiterais tout d'abord remercier mes directeurs de these. Un grand merci a
Samuel, qui il y a 4 ans maintenant a accepté de travailler avec un presque economiste
presque statiscien (moi-même du coup) pour faire du machine learning avec des
modèles de climat. Le pari était risqué et je l'en remercie le plus sincèrement du
monde de l'avoir pris. Pendant ces 4 années il a eu la patience de me former aux
enjeux de la modélisation climatique, de me transmettre un bout de son experience
sur les modèles regionaux et la curiosité de se confronter aux obscures réseaux de
neurones. Je le remercie aussi d'avoir pris le temps comprendre ce qu'il se passait
dans ma tête (ou du moins de faire semblant) et de m'avoir fait confiance pour
qu'au final on arrive à quelque chose de pas trop mal. Enfin merci d'avoir su gerer
les différentes solliciations parce que je n'aurais pas su le faire et cette thèse n'aurait
jamais fini.

Merci aussi à Sébastien qui a accepté de diriger cette thèse alors que le sujet était encore plus loin de son champ d'expertise. Au delà de son soutient et ses conseils au combien importants, il se sera attaché durant ces 3 années à me transmettre la rigueur mathématique à laquelle l'impatient que je suis aura eu du mal à se plier.

Dans l'ordre de ceux qui m'ont mis sur ce chemin périlleux de la thèse et de la
recherche climatique je remercie Bastien qui m'a montré qu'après la TSE une autre
vie était possible. Je tiens à remercier tres chaleureusement Hervé qui m'a ouvert
les portes du CNRM au printemps 2018 et qui a poussé pour faire de moi le premier
alternant du labo. Vient ensuite Aurélien, qui après 9 mois de collabration que je
trouvais interessante a décidé d'aller encadrer et former ailleurs. Je le remercie,
quand même, d'avoir ardé un œil sur nos avancées et du temps pour des discussions
et conseils toujours très utiles.

Je tiens à remercier tout paticulièrement Antoinette, Lola, Cécile, Pierre, Flo-
rence qui ont contribué à faciliter cette thèse en étant toujours disponibles et d'une
aide précieuse face à mes récurentes (et souvant redondantes) questions techniques.
J'en profite ici pour remercier tou.te.s les membres de la MOSCA-team pour leur
acceuil et leur bonne humeur permanente.

Plus généralement je tiens à remercier l'ensemble des collègues grâce à qui il est
plutôt agréable de venir au labo. Ici je pourrai faire une longue liste de noms en
esperant n'oublier personne mais je vais en fait essayer d'être efficace parce que c'est
suffisamment long comme ca et dire en une phrase (longue certes) un grand merci
à tout GMGEC pour ces moments partagés à la pause café (j'en profite ici pour
addresser un merci particulier à tou.te.s les faiseuses et faiseurs de café, Christine et
Florence en tête), autour d'un délicieux repas à Eurest ou encore au bar à partager
de la boisson alcoolisé (et là je remercierai spécialement les plus fidèles car ils ont
largement contribué au bon déroulé de cette thèse: Robin, Erwan, Raphael, Nicolas,
Maya, William, Romain, Nils, Chloé, Gildas, Thomas...).

Enfin, last but not least, je remercie tous mes proches qui ont été un soutient
indispensable à cette thèse. Mes parents qui m'ont toujours permis d'avancer dans
les meilleurs conditions et toujours soutenus dans mes différents projets. Mes frères
qui ont bien de la chance de m'avoir comme frère. Mes amis, qui ne liront proba-
blement jamais ces lignes et ils ont bien raison, mais heureusement qu'ils étaient là
pour me sortir la tête du guidon. Et je terminerai par celle que cette thèse a le plus
impacté car elle m'a supporté au quotidien pendant 3 ans, et je suis bien placé pour
savoir que (notamment ces derniers mois) c'était pas facile, un très grand merci à
Anja d'être toujours et pour longtemps encore inshallah à mes cotés.
# Contents

## 1 Introduction 10

1.1 What is Climate? 10

1.1.1 Definition 10

1.1.2 Climate System 11

1.1.3 Natural variability and climate change 12

1.2 Climate modelling 14

1.2.1 What is a climate model? 14

1.2.2 Climate simulations 15

1.2.3 Different sources of uncertainties in climate projections. 16

1.3 Regional climate change information 20

1.3.1 Dynamical downscaling 20

1.3.2 Statistical Downscaling 30

1.4 Objectives and scientific questions of this study. 35

## 2 Why deep neural networks are a credible strategy to build the RCM-emulator ? 38

2.1 A neural network 38

2.2 Training the network: the backpropagation algorithm 40

2.3 The last milestones: the ReLu function and the computational resources 44

2.4 Convolutional neural network (CNN) 44

2.5 Other great families of architecture 48

2.6 Deep Learning and climate applications 49

2.7 Link with this work 52

## 3 Regional Climate Model emulator based on deep learning: concept and first evaluation of a novel hybrid downscaling approach 53

3.1 Introduction 54

3.2 Methodology 56

3.2.1 Models and simulations 56

3.2.2 Predictors 57

3.2.3 Predictands 59

3.2.4 Deep learning with UNet 60

3.2.5 Training of the emulator : Perfect Model Framework 63

3.2.6 Evaluation metrics 64

3.2.7 Benchmark 65

3.3 Results 66

3.3.1 Computational efficiency 66

3.3.2 Evaluation step 1 : Perfect model world 66
<table>
  <tbody>
    <tr>
      <td>3.3.3</td>
      <td>Evaluation step 2: GCM world</td>
      <td>73</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Discussion</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>On the inconsistencies between GCM and RCM</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.4.2</td>
      <td>On the stationary assumption</td>
      <td>79</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>On the selection of the predictors</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.4.4</td>
      <td>On the non reproducibility of neural network training</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Conclusion</td>
      <td>83</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Appendices</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.6.1</td>
      <td>Emulators inter-comparison</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.6.2</td>
      <td>Description of the UNet-emulator architecture</td>
      <td>89</td>
    </tr>
    <tr>
      <td>3.6.3</td>
      <td>Supplementary Material: Summary tables</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Study of the RCM-Emulator applicability to GCM ensembles</td>
      <td>92</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Introduction</td>
      <td>93</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Methodology</td>
      <td>95</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>The ALADIN63 matrix</td>
      <td>95</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Recall on the main emulator and Experiment design</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Justification of the perfect model framework</td>
      <td>98</td>
    </tr>
    <tr>
      <td>4.2.4</td>
      <td>Metrics, Benchmark and results presentation</td>
      <td>99</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Perfect Model evaluation</td>
      <td>101</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Main emulator</td>
      <td>101</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Cross-Validation on other training sets</td>
      <td>111</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Application on GCM data</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Discussion</td>
      <td>120</td>
    </tr>
    <tr>
      <td>4.5.1</td>
      <td>On the standardization method</td>
      <td>120</td>
    </tr>
    <tr>
      <td>4.5.2</td>
      <td>Some tests on inputs choices</td>
      <td>121</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Conclusion</td>
      <td>122</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Precipitation downscaling with the RCM-Emulator: The role of the loss function</td>
      <td>124</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Introduction</td>
      <td>125</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Methodology</td>
      <td>127</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>Datas: the ALADIN63's matrix</td>
      <td>128</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>RCM-Emulator conception: Predictors, architecture and target</td>
      <td>128</td>
    </tr>
    <tr>
      <td>5.2.3</td>
      <td>Loss function for the neural network training</td>
      <td>130</td>
    </tr>
    <tr>
      <td>5.2.4</td>
      <td>Evaluation Metrics</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Perfect model Evaluation</td>
      <td>138</td>
    </tr>
    <tr>
      <td>5.3.1</td>
      <td>Comparison of the three emulators</td>
      <td>138</td>
    </tr>
    <tr>
      <td>5.3.2</td>
      <td>Deeper evaluation of Emul-ASYM</td>
      <td>146</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>GCM data application</td>
      <td>155</td>
    </tr>
    <tr>
      <td>5.4.1</td>
      <td>Illustration of the daily GCM/RCM differences</td>
      <td>155</td>
    </tr>
    <tr>
      <td>5.4.2</td>
      <td>Present climate analysis</td>
      <td>157</td>
    </tr>
    <tr>
      <td>5.4.3</td>
      <td>Climate change analysis</td>
      <td>157</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Discussion</td>
      <td>160</td>
    </tr>
    <tr>
      <td>5.5.1</td>
      <td>On the asymmetric loss</td>
      <td>160</td>
    </tr>
    <tr>
      <td>5.5.2</td>
      <td>On the computational efficiency</td>
      <td>162</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Conclusion</td>
      <td>163</td>
    </tr>
  </tbody>
</table>
<table><tr><td>6 Conclusion & Perspectives</td><td>166</td></tr><tr><td>6.1 RCM-Emulator concept</td><td>167</td></tr><tr><td>6.2 Efficiency</td><td>168</td></tr><tr><td>6.3 Trustability</td><td>169</td></tr><tr><td>6.4 Applicability step 1: different sources of input</td><td>171</td></tr><tr><td>6.5 Applicability step 2: GCM world</td><td>172</td></tr><tr><td>6.6 More general perspectives</td><td>173</td></tr></table>
# Acronyms

<table><tr><td>IPCC</td><td>Intergovernmental Panel on Climate Change</td></tr><tr><td>CMIP</td><td>Coupled Model Intercomparison Project, different phases exist, the last is CMIP6</td></tr><tr><td>AR</td><td>Assessment Report, they are the IPCC reports 6 have published by now.</td></tr><tr><td>GCM</td><td>Global Climate (or Circulation) Model</td></tr><tr><td>RCM</td><td>Regional Climate Model</td></tr><tr><td>CPM</td><td>Convection Permitting Model</td></tr><tr><td>RCP</td><td>Representative Concentration Pathway (CMIP5 scenarios of greenhouse gas emissions)</td></tr><tr><td>SSP</td><td>Shared Socio-economic Pathway (CMIP6 scenarios of greenhouse gas emissions)</td></tr><tr><td>CORDEX</td><td>Coordinated Regional climate Downscaling Experiment</td></tr><tr><td>DD</td><td>Dynamical Downscaling</td></tr><tr><td>ESD</td><td>Empirical Statistical Downscaling</td></tr><tr><td>PP</td><td>Prefect Prognosis</td></tr><tr><td>MOS</td><td>Model Output Statistic</td></tr><tr><td>WG</td><td>Weather Generator</td></tr><tr><td>GHG</td><td>GreenHouse Gas</td></tr><tr><td>TAS</td><td>Surface Air Temperature</td></tr><tr><td>PR</td><td>Precipitations</td></tr><tr><td>UA,VA</td><td>Wind vector components</td></tr><tr><td>HUS</td><td>Specific Humidity (sometimes Q)</td></tr><tr><td>Z</td><td>Geopotential altitude</td></tr><tr><td>SST</td><td>Sea Surface Temperature</td></tr><tr><td>NN</td><td>Neural Network</td></tr><tr><td>CNN</td><td>Convolutional Neural Network</td></tr><tr><td>RNN</td><td>Recurrent Neural Network</td></tr><tr><td>LSTM</td><td>Long Short-Term Memory</td></tr><tr><td>GAN</td><td>Generative Adversarial Network</td></tr></table>
# Chapter 1

# Introduction

## 1.1 What is Climate?

### 1.1.1 Definition

*Climate is what you expect, weather is what you get.*
*Andrew John Herbertson*

Defining climate can be delicate. The best way to understand it is maybe to differentiate it from the weather. Climatology and meteorology are two different disciplines looking at the same object but at different time scales. The object is the Earth system and, more specifically, its atmosphere component. Meteorology uses the current knowledge of the atmospheric processes and the observation of the previous and actual atmospheric conditions to eventually predict the weather for the next hours and days. The system is by nature chaotic (Lorenz, 1963). The smallest differences in the *initial conditions*<fnref n="1" /> can lead to a radically different evolution of the atmospheric circulation. The chaotic nature of the atmospheric system limits the weather prediction to a maximum of two weeks.

The climatology goes beyond these short time variations. Climatology is directly linked to statistics and probabilities. It aims to comprehensively describe the distribution of the possible weather that can be encountered in a given region and period of time. It does not look at the weather of a given day but at all weathers observed over decades in a region to understand the mean state of the local climate, its variability, its seasonal evolution or the extremes of climate events, for example. More practically, the climate, such as the weather, is described through variables and quantities like temperature, humidity, wind or cumulative precipitations. The typical period length to study the climate is 30 years (as defined by the World Meteorological Organization). It allows having enough samples of weather to obtain some robust statistical properties. The climate of a region depends on many factors like the latitude, surrounding topography, distance to the sea, presence of a lake or forest etc..

In this context, climate change is the transformation of the statistical distribution of the variables describing the climate. Proving it is establishing the statistical

<footnote n="1">¹i.e. the atmospheric conditions at a given time *t* marking the beginning of the experiment.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.1:** Climate system representation, extracted from IPCC 4th Assessment Report<fnref n="2" /> (AR4 Solomon et al, 2007).</figcaption></figure>

difference between two periods of time, taking into account the internal variability of the climate due to its chaotic nature.

## 1.1.2 Climate System

The climate system is the interactive system of 5 major components: the atmosphere, the hydrosphere, the cryosphere, the lithosphere and the biosphere. This whole system is directly dependent on external forcings<fnref n="3" /> like the sun, which is the most influential. It is the primary source of energy on Earth. Because it affects the planet unequally, it creates temperature gradients, putting the atmosphere and the ocean in movement to redistribute this heat. The poles are cooler because they receive less radiation but also because the ice reflects the radiation from the sun. This significant gradient between the poles and the tropics creates atmospheric circulation. The land surfaces (lithosphere) also play an essential role as they evaporate water differently, absorb or reflect solar radiation, or because of the topography, all these interactions impact atmospheric circulation.

<footnote n="2">²IPCC Assessment Reports are the well known reports about knowledge on climate change and its impacts. IPCC is an intergovernmental organisation to which thousands of people contribute. The volunteer authors of those reports review the thousands of scientific papers published each year to provide a comprehensive summary of what is known about the drivers of climate change, its impacts and future risks, and how adaptation and mitigation can reduce those risks.</footnote>

<footnote n="3">³External agent that impacts the system while being outside.</footnote>
The global ocean covers 71% of the Earth. It is the primary climate regulator through permanent interactions with the atmosphere. It receives, stocks and transports the sun's energy. The oceanic circulation transports heat from the tropical regions to the higher latitudes. The cold water from these regions gets denser and dives into the deep ocean. Those mechanisms combined with the surface wind are the start of the reasons behind ocean circulation. The time scale of these transfers varies from season or year for the tropical regions to centuries for the profound water changes. Oceans also play a significant role in the water cycle. Heated by the sun, the surface water evaporates continuously, is transported by atmospheric currents, condensates in clouds, and precipitates over land or ocean. Freshwater exchanges with the atmospheres, land surfaces or the cryosphere create a gradient of density in the ocean, also contributing to the ocean circulation. Because of its much slower circulation, the ocean generally impacts the atmospheric circulation at longer times scales.

The atmosphere contains a mix of gases. Among them, some *trace gases* interact with the sun's radiations. This is the greenhouse effect. The atmosphere directly reflects a third of the sun's radiations in space. The system absorbs the rest, which is then re-emitted by the atmosphere and the land surface. The greenhouse gases (GHG, e.g. carbon dioxide, methane, ozone or water vapour) absorb the infrared radiation from inside and outside the system and re-emit them in all directions. These radiations are then trapped in the system and warm the planet to a temperature adapted to support life.

### 1.1.3 Natural variability and climate change

Since the work of Lorenz (1963), it is common knowledge that the atmosphere is a chaotic system. It is sensitive to any change in the initial conditions making it the most unstable component of the climate system. The climate research community often refers to this chaotic aspect of the system as the natural or internal variability of the system in opposition to the external forcing.

Natural climate variability encompasses various processes driving the variability at different time scales. Williams et al (2017) retraces the different kinds of mechanisms behind internal atmospheric variability. It goes from the turbulence induced by the diurnal cycle or rapid air-sea interactions, which have immediate impacts on the atmosphere, to multiyear oscillations coupling atmosphere and ocean circulations such as the El-Nino phenomena. The critical point here is that naturally, the climate system greatly varies, following different modes impacting local climates differently and at different time scales. This chaotic nature has strong impact on the way to predict weather and to model climate requiring probabilist or ensembles-based approach.

On top of the atmospheric variability, some natural external forcing might strongly influence the system. Volcanoes' eruptions emit gases absorbing the sun's radiations and cooling the atmosphere. For example, The Mount Pinatubo eruption in 1991 caused a cooling of the global temperature of half a degree for the following months. The variations in the sun radiations received by the Earth due to the sun itself or
changes in the Earth’s orbit also impact climate as they modify the quantity of energy received by the system.

Human activity also influences the climate system and plays the role of anthro-
pogenic external forcing. Modifying land surfaces for agriculture or creating large
urban areas has apparent consequences on the local climate but can also have reper-
cussions at larger spatial scales. Since the industrial revolution, human activities
have emitted a considerable amount of greenhouse gases. The concentration of car-
bon dioxide in the atmosphere increased by 150% compared to the levels of 1750,
methane has more than doubled (+ 260 %) and the nitrous oxide increased by 123%
(cf. Figure 1.2). The increase in greenhouse gas concentration modifies the compo-
sition of the atmosphere. It warms the planet as it reinforces the greenhouse effect<fnref n="4" />.
The speed of the atmospheric changes induced by human activities is unprecedented
in Earth’s history and has alarmed the scientific community for decades. The inter-
national community decided to merge efforts to evaluate the impact of this GHG
increase on our climate. The main tools used to do so are climate models.

<figure><img src="image_4.png" /><figcaption>**Figure 1.2:** Evolution of the main greenhouse gas concentration over million year. IPCC AR6, Masson-Delmotte et al (2021), chap5. The reader should notice here that the x-axis is not linear.</figcaption></figure>

<footnote n="4">⁴The warming process due to GHG emissions is more complex, in particular including positive feedback loops, but we keep it simple here as it is not the purpose of this manuscript.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Scheme of a GCM.</figcaption></figure>

## 1.2 Climate modelling

Climate models are numerical tools representing the climate system. The interest in such models is double. They allow a better understanding of the system by realising numerical experiments impossible to run in the real world. Furthermore, they permit the exploration of other climates by modifying some components to reconstruct past periods or to make projections about future climates.

### 1.2.1 What is a climate model?

There is a wide variety of climate models. Simple climate models represent part of the system at the global scale, like the system's energy balance. Such models are used to assess the global temperature according to the CO2 concentration in the atmosphere, for example. In this work, I will focus only on more complex models that aim to represent the entire system. The starting idea of such models is to represent the time evolution of the atmosphere.

The atmospheric circulation is represented by the equations of fluid mechanic theory, respecting the physical laws of mass and energy conservation. The atmosphere is divided into small cubes (or meshes or grid points) in which the equations are solved from one time step to the other as well as interactions with the surrounding cubes. The grid is three-dimensional in order to represent both horizontal and vertical exchanges. At each time step and in each mesh, the model returns climatic variables like temperature, humidity or wind. A similar numerical representation of
oceanic circulation is done with a dedicated model. The atmospheric and oceanic
models are then coupled to represent the interactions between the two components.

Following the same principle, land surfaces, continental hydrology cycle and sea
ice exchanges with the atmosphere and ocean are represented by dedicated mod-
els and coupled with the atmosphere and ocean components. These are the main
components of the so-called General Circulation Models or Global Climate Models
(GCMs). Besides, the model receives external forcing, like the truth climate sys-
tem. The sun radiations entering the system or the greenhouse gas concentrations
or emissions resulting from anthropogenic activities are given as input to the model.

The complexity of GCMs is constantly increasing to represent better each phys-
ical or chemical process playing in each system component. The most recent and
complex models are called Earth System Models. They target to represent the entire
carbon cycle. For example, they represent the life of aerosols (suspended particles
in the atmosphere): their formation, transportation by wind, and interaction with
radiations, clouds or precipitation. The representation of land and biosphere is also
more accurate: different types of vegetation are represented as they play a signifi-
cant role in the system through photosynthesis or evapotranspiration processes.

The cube size is a determinant parameter in a climate model. As mentioned earlier, some climate processes can be highly local such as the formation of clouds or the occurrence of rain. Therefore, the smaller the cube, the finer is the representation of the system since the fluids mechanic equations are more precisely resolved. However, small cubes imply more cubes, so more equations to solve, more variables to estimate and save, and more computational resources. The model resolution (the grid point size) results from a compromise between the precision of the model, its complexity<fnref n="5" /> and the available computing resources. GCMs have mostly a horizontal resolution of about one to two hundred kilometres. It implies that some mechanisms cannot be simulated. *Empirical* or *physical parametrisations* are set to represent the effects of not represented processes on the model variable at the meshes scale. The resolution of the model grid and the choices of those parameters explain a large part of the differences between climate models.

1.2.2 Climate simulations

Climate models can be used for different purposes. If we put aside a better
understanding of the processes, there are two main applications of climate models:
exploring the past or the future. First, it is used to create climate *reanalysis*. They
reconstruct the recent climate history by combining climate models and all available
observations (in-situ records and satellite). They the best estimate of the true
atmospheric circulation of the last decades. They are widely used in the scientific
community since they are continuous and homogeneous in space and time and match
local observations. They are notably used to evaluate climate models, especially
for non-observed variables. In a further extent, Global Climate Models are used
to reconstruct past climates for which no records exist. This is paleoclimatology,

<footnote n="5">⁵i.e. the number of climate system components that are represented.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 1.4:** Illustration of two types of models ensembles exploring two different sources of uncertainty.</figcaption></figure>

combining climate models and the study of geological sediments. Studying ancient climates is necessary to understand the situation where we are know.

Second, climate models are widely used to create future climate projections. The idea is to study the effects of human activities and related evolving emissions on the system. Such simulations generally start at the beginning of the industrial revolution (1850). They are pushed to a few centuries after today, even if the most studied period is the end of the current century. The study of climate change compares a future climate with either pre-industrial or more recent climates (generally 1981-2010). Nevertheless, there is considerable uncertainty from multiple sources associated with such simulations.

### 1.2.3 Different sources of uncertainties in climate projections.

Climate models are far from being perfect. Even in the most complex model, it is impossible to represent all processes, bringing some uncertainty to the accuracy of the representation. We also mentioned above the role of resolution or parameters choices in the differences between models. The answer to this first source of uncertainty is the multiplication of the number of models. Many global climate models are developed in many research centres worldwide. Each of them is wrong in its way. Since the mid-'90s, the scientific community decided to merge the efforts and build
<figure><img src="image_7.png" /><figcaption>**Figure 1.5:** Time evolution of the main greenhouse gas concentration in the 4 RCP scenarios. From AR5, Tignor (2013).</figcaption></figure>

great models intercomparison exercises called CMIP (Coupled Model Intercomparison Project). The aim is to build large ensembles of climate simulations from the existing models. Comparing models built by different teams who made different choices brings robustness. It helps to identify the weaknesses and the strengths of each model. Moreover, if all models agree on a given message, it gets more trustable. Figure 1.4 illustrates the time series of a multi-model ensemble. In Figures 1.6 and 1.7, commented later on, the straight line shows the multi-model mean and the shaded area represents a measure of the model ensemble spread.

The second source of uncertainty is brought by the internal (or natural) variability of the climate and so of the models. Indeed, climate models aim to reproduce the chaotic dynamic of the climate system, so the simulations are highly dependent on the initial conditions fixed at the beginning of the simulation. In order to take into account the possible impact of the internal variability, the same climate models are run multiple times with the same configuration and external forcing but different initial conditions. These multiple simulations are members of the model ensemble. To a larger extent, SMILEs (Single Model Initial-condition Large Ensemble) are vast ensembles of climate model members (more than 30). The number of studies using such ensembles has primarily increased over the last decade, and they help to quantify the model's internal variability. Maher et al (2021) proposes an overview of the application of SMILEs.
<figure><img src="image_8.png" /><figcaption>**Figure 1.6:** Global surface temperature evolution in the CMIP5 ensemble. From Knutti and Sedláček (2013).</figcaption></figure>

The last source of uncertainty concerns the future evolution of human activities. In order to answer this uncertainty, the climate community established different scenarios of greenhouse gas emissions. The scenarios are determined according to key drivers such as demography, economy trajectories etc. Following these scenarios, tables of greenhouse gas concentrations are used to drive the climate models.

CMIP exercises comport several members of several GCMs driven by different scenarios of human emissions. CMIP5⁶ is the fifth generation of global climate models ensemble used for the fifth report of the IPCC published in 2014, Tignor (2013). CMIP5 regroups about 50 climate models, resulting into about 1000 long-term simulations (i.e. at least 20th and 21st centuries). Four main scenarios were chosen to drive CMIP5 models after 2006, while the historical part is driven by a reconstruction of historical greenhouse gas emissions. The RCP2.6, RCP4.5, RCP6.0 and RCP8.5 scenarios have been introduced in Meinshausen et al (2011) and are named accordingly to the level of global radiative forcing induced by the GHG concentrations they reach in 2100. The four scenarios have similar trajectories until 2030 and then diverge, as shown in Figure 1.5. The RCP2.6 is the most “optimistic” one, as it peaks around 2030 and then projects a slow reduction of the CO2 concentration while the methane concentration diminishes drastically to reach mid-20th century levels before 2100. The RCP8.5 corresponds to an exponential increase in GHG until the end of the 21st century. The two others are intermediate scenarios. Figure

<footnote n="6">⁶The models used in this work are from CMIP, which is why CMIP5 scenarios and figures are presented here.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 1.7:** Illustration of the time evolution of the global surface temperature in CMIP6 model with and without anthropogenic greenhouse gas emissions. From IPCC AR6, Masson-Delmotte et al (2021).</figcaption></figure>

1.6 shows the global surface warming projected for the different scenarios. The solid line shows the model's mean, and the shaded area shows the model spread. CMIP6 is the most recent CMIP exercise and had been used in the sixth IPCC assessment report (AR6, 2021/22, Masson-Delmotte et al, 2021). It regroups about 70 models and more than a thousand simulations. New scenarios have been defined under “Shared Socioeconomic Pathways” (SSPs).

Besides historical and scenario-based projections, CMIP exercises include other kinds of experiments. For example, one consists in running the model on the historical period (1850-2006 for CMIP5, -2014 for CMIP6), removing the GHG emissions due to human activities. Figure 1.7, extracted from the AR6 summary for policy-makers, shows the experiment's results on the CMIP6 models. This figure is a good illustration of the use of large ensembles. Here all models agree that the observed global temperature warming is not due to the climate's internal variability or natural forcings. The observed temperature time series fits in the range of simulations as a single realisation among the possible ones. The variety of models ensembles establishes the range of possibilities.

The repartition of the total uncertainty among the three sources we defined before is an essential topic with multiple studies trying to establish it. Lehner et al (2020), following the framework set in Hawkins and Sutton (2009), assessed this question using seven SMILEs models and both CMIP5 and six ensembles. Figure 1.8 shows the study's results for the global scale temperature. Over the first decades, internal variability is the primary source of uncertainty. It is followed by the model
uncertainty and scenario uncertainty, which becomes the leading one from the mid-
century. We can notably see that internal variability plays no role in the global
temperature uncertainty after some decades.

**1.3 Regional climate change information**

*Various mechanisms operating at different time scales can modify the amplitude of the regional-scale response of temperature, and both the amplitude and sign of the response of precipitation, to human influence (high confidence).*

IPCC AR6, Working Group I, (Masson-Delmotte et al, 2021), Chapter 10.

If climate change is a global phenomenon, its impacts vary from one region to the other. Worldwide, the consequences of global warming are already impacting societies. In order to conduct studies on the impacts on agriculture or hydrological cycle, for example, it is crucial to produce climate change information at regional scales. The regional scale varies from a few thousand kilometres for sub-continental regions such as the Mediterranean basin to some kilometres representing cities, valleys or coastal climates. In the rest of this manuscript, we will also talk about local scales to refer to those small areas that share the same climate and where climate change information is needed to set up adaptation strategies.

Global climate models have a resolution of about 100 to 200 km. Several studies have demonstrated the difference between nominal resolution and effective resolution. The nominal resolution is the one we defined earlier; it corresponds to the actual size of the grid points. The effective resolution of the model is the resolution where the model is reliable. Klaver et al (2020) concludes that the effective resolution of a model is about 3 to 5 times its nominal resolution. This is too coarse to account for local drivers of regional climate correctly. Thus, GCMs provide, at best subcontinental information on climate change. For example, in the last CMIP exercises, models agree on more precipitation over Northern Europe in winter and a drier summer over southern Europe at the end of the century for the RCP4.5 scenario and above.

Nevertheless, impact studies and adaptation policies necessitate much higher reso-
lution. The higher resolution brings a better representation of processes parametrised
in GCMs, and a better representation of local drivers such as mountains, coastlines,
towns, great lakes or islands. The solution brought by the community is to downscale
the GCM simulations. Two great families of downscaling have appeared: dynamical
and statistical downscaling.

**1.3.1 Dynamical downscaling**

Dynamical downscaling relies on another kind of climate model called Regional Climate Models (RCMs), which focus only on a regional domain. Regional climate modelling is quite a recent activity as it appeared about 30 years ago. Dickinson et al (1989) and Giorgi and Bates (1989) introduced the first RCM in the late 80s. Since
<figure><img src="image_10.png" /><figcaption>**Figure 1.8:** Figure and legend from Lehner et al. (2020): (a–c) 10-year running means of global annual mean temperature time series from (a) SMILEs, (b) CMIP5 and (c) CMIP6, with observations (Rohde et al., 2013) superimposed in black, all relative to 1995–2014. For SMILEs, the ensemble mean of each model and the multimodel average of those ensemble means are shown; for CMIP the polynomial fit for each model and the multimodel average of those fits are shown. (d–f) Illustration of the sources of uncertainty in the multimodel multi-scenario mean projection. (g–i) Fractional contribution of individual sources to total uncertainty. Scenario uncertainty in SMILEs in (g) is taken from CMIP5, since not all SMILEs offer simulations with multiple scenarios. (d–i) In all cases, the respective multimodel mean estimate of internal variability (Imean) is used.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 1.9:** Scheme of dynamical downscaling process, from Giorgi (2019)</figcaption></figure>

then, the number of studies including an RCM never stopped increasing. Hereafter, we review some key conceptual and technical points of RCMs.

## Concept and technique

Regional climate models are numerical models based on the same concept as global models, except that they run on a limited area domain such as Europe (cf Figure 1.9)<fnref n="7" />. It allows them to run at higher resolution, between 10 and 50 km. The regional climate model is embedded or nested in a GCM simulation to ensure the transfer of information with the global scale. We say that the GCM drives or forces the RCM. Figure 1.9 illustrates the nesting scheme. RCMs generally have a temporal resolution of less than half an hour. At each time step, the RCM receives information about general circulation outside its domain from the GCM simulation thanks to a set of lateral boundary conditions (LBC). LBCs are contained in the main atmospheric variables: the vertical distribution of moisture, temperature, zonal and meridional wind components in the vertical, and surface pressure. They are extracted from the GCM driving simulation and transmitted to the RCM in the so-called buffer zone (cf Fig. 1.9) thanks to some relaxation technique (Davies and Turner, 1977). LBCs are generally stored every 6hrs in the GCM simulation and interpolated to match the RCM temporal and spatial resolution. The sea surface temperature (SST), the sea ice cover, the aerosols and greenhouse gas forcings from the driving GCM simulation are also given as input to the RCM.

RCMs have a higher resolution than the driving GCM. The first RCM was introduced to downscale a 500km resolution GCM over the Yucca mountains in the

<footnote n="7">⁷An other family of RCMs are global models with a stretched grid allowing a higher resolution over a specific domain. However, for simplicity we do not consider them here</footnote>
United States of America in the late 90s. It had a horizontal resolution of 70km. Nowadays, thanks to a significant improvement in computational efficiency, most GCMs have a resolution of about 100 to 200km, with some reaching higher resolutions close the 50km, with computational costs still allowing to run numerous long-term simulations. RCMs simultaneously increased their resolution to 10 to 25km. The ratio between the GCM and RCM resolution is an important parameter to ensure the accuracy of the downscaling, and it should stay below 12 times (Denis et al, 2003). Moreover, RCMs include their land surface scheme at a higher resolution than the GCM to take better account of local land surface characteristics (reliefs, littoral, land-use).

The nesting technique of the RCM is a widely debated topic, and various methods have been developed. The one described above is called one-way nesting. The RCM receives at its boundary regular information from the GCM but does not feed back to it. Inside its domain, the RCM develops its internal physics and dynamics. However, at large scales<fnref n="8" /> the RCM should not deviate from the GCM simulation. It is designed to add realistic sub-GCM grid scale details but not to correct or modify the GCM global circulation. It implies using good GCMs in order to get good high-resolution RCM simulations. However, it depends strongly on the domain size. In a large domain, the boundary control is less strong, and the RCM is more susceptible to developing its own large scales (Laprise et al, 2008; Miguez-Macho et al, 2004).

A variant of the classical one-way nesting allowing controlling the eventual RCM deviation is called spectral nudging (Waldron et al, 1996; Von Storch et al, 2000). Such a method constrains the RCM stronger, ensuring a better consistency with the GCM large scale. It is the same as one-way nesting but also controls the large scales on top of the domain. Above a given altitude, the RCM atmosphere large scales are controlled to follow the GCM information for some key variables (e.g. temperature, moisture or wind components). The soundness of spectral nudging has been widely debated in the literature. Some studies argue that it limits the RCM internal processes too much, creating inconsistent dynamics or failing to represent some local extremes (Alexandru et al, 2009; Separovic et al, 2012; Rummukainen, 2010; Castro, 2005). Even if multiple studies have shown that with the appropriate settings, the spectral nudging technique increases the consistency with the driving GCM simulation without evident deterioration of the local processes (Sanchez-Gomez et al, 2009; Lucas-Picher et al, 2016; Colin et al, 2010), most RCM simulations are performed following the classical one-way nesting methodology. Two-way nesting simulations have also been performed with the RCM feeding back to the GCM simulation. Such experiments did not show a strong regional to global impact, except for specific regions like the poles (Giorgi, 2019) and the tropics (Lorenz and Jacob, 2005). Besides, their extremely high costs are a severe brake on their use.

Such as GCMs, the development of RCM goes in three directions: a better

<footnote n="8">⁸The *largest scales* are challenging to define. We consider them in this manuscript to match the effective resolution of the GCM, where the GCM is doing a good job. However, the closer we go to the surface, the more the resolution impacts the quality of the representation. So we consider the large scales in opposition to the local scales where the RCM aims to add value. When we refer to large scales in this manuscript, we mean the high altitude general circulation of the atmosphere represented by variables at low resolution and low-pressure levels.</footnote>
model content (improved parametrisation), a better consideration of the different
components of the climate system and an increase in resolution to represent better
the most local processes. There has been a strong development of Atmosphere-Ocean
coupled RCM with some including biogeochemistry representation in both land and
ocean (called RCSM or RESM, Somot et al, 2008; Sevault et al, 2014; Sitz et al,
2017; Sein et al, 2015). Interactive models for the representation of aerosols have
also been developed as these particles have a substantial impact on local climates
and global warming (Nabat et al, 2014). On the other hand, the last decade has
seen the apparition of Convection Permitting Models (CPMs, or CP-RCMs) (Prein
et al, 2015; Ban et al, 2021; Lucas-Picher et al, 2021; Caillaud et al, 2021). They are
regional climate models running at a resolution between 2 and 5km, and allowing
the numerical solving of the deep convection, responsible for local storms. In order
to respect a reasonable resolution ratio, CPMs are most of the time nested into
RCMs simulations. It is then a two-step nesting.

RCMs ensembles

Regional Climate Models constitute a novel source of uncertainty. They are a
modelling step more with different choices in the nesting method, parameters set-
tings or the complexity of the model. Such as GCMs, the choice of the RCM is then
a source of uncertainty. In order to take into account all sources of uncertainties, it
is necessary to create RCMs ensembles or matrices, where several GCMs simulations
drive several RCMs over the same regional domain and period (Déqué et al, 2012).
Because of the GCM uncertainties, it is necessary to take GCM simulations driven
by different scenarios and several members of the same Scenario-GCM pairs. It
leads to 4 dimensional [Scenario, GCMs, Member, RCM] matrices. Moreover, inter-
nal variability plays a preponderant role at regional scales and regarding some local
and extreme events. Ideally, hundreds to thousands of RCMs simulations would be
needed in a given region to cover those uncertainties at the maximum and therefore
provide a reliable local climate information. However, here comes the most signifi-
cant limitation of dynamical downscaling. Even if they have a limited area domain,
the high resolution implies much more grid points and much higher computational
costs. For example, a global grid at 150km resolution includes about 30 thousand
points. In comparison, a large European domain at 12km includes about 200 thou-
sand of them, which must be multiplied by the altitude levels (between 50 to 100).
Moreover, RCM have a higher temporal resolution which also increase their costs.
The number of long-term RCM simulations is much more limited than the GCMs
with the existing computational resources.

Facing this issue, the scientific community, especially in Europe, created ensem-
ble projects similar in spirit to the CMIP exercises. PRUDENCE (Déqué et al,
2007) and ENSEMBLES (Déqué et al, 2012) projects successively contributed to
the creation of the first European RCM-GCM matrices at 50km and 25km resolu-
tion. Similar exercises have been raised for the different continents leading to the
CORDEX initiative (Giorgi and Gutowski, 2015). It aims to coordinate the Re-
gional Downscaling activities worldwide to provide a coherent and reliable source of
downscaled information. The different regional communities regroup around some
well-defined domains. For each domain, the participants agree on a simulation pro-
tocol to coordinate the production of downscaled simulations. For example, GCMs
<figure><img src="image_12.png" /><figcaption>Figure 1.10: EURO-CORDEX matrix, October 2022.</figcaption></figure>

are selected according to their performance and spread of future projections. Each participant picks some GCMs to downscale among this list. The idea is to explore at maximum among the available and realistic GCMs. A minimum number of scenarios to downscale or the allowed nesting technique are also defined in advance. The aim of such exercises is, for a given region, to coordinate the filling of RCM-GCM-Scenario matrices to explore the possible evolution of local climates.

The production of CORDEX simulations is done in three steps. First, all participating RCMs perform a perfect condition simulation with LBCs from reanalysis over a recent period (1980-present). This simulation, called evaluation simulation, is essential to evaluate the RCMs against observations. Second, the RCMs are used to downscale the historical (1950-present) simulations from the selected GCMs. It helps validate the RCM driven by GCM simulations by comparing the climatology and variability with the evaluation run. It also serves as a reference to assess the eventual changes induced by global warming, as the third steps consist of downscaling GCM projections under various scenarios over a future period running until 2100.

EURO-CORDEX is the most collaborative community and filled up the biggest matrix (cf Fig 1.10 Vautard et al, 2020; Coppola et al, 2021). It includes 13 RCMs, 8 GCMs from CMIP5 and 3 RCP scenarios (2.6, 4.5 and 8.5). After more than ten years of coordinated efforts on the potential 312 projection simulations, only a hundred have been performed. On this number, we should also account for the historical and evaluation simulations, which lead to a number of about 170 simulations. This is largely more than any other CORDEX domain, but still, it is insufficient. Indeed the matrix is highly unbalanced, and only a fifth of the existing CMIP5 GCMs are considered. Evin et al (2019) showed that an unbalanced or too weak matrix could lead to a wrong estimation of future changes.

Moreover, in these exercises, internal variability is rarely explored. In the EURO-
<figure><img src="image_13.png" /><figcaption>**Figure 1.11:** Sub-continental uncertainty decomposition in RCM ensemble, from Evin et al (2021).</figcaption></figure>

CORDEX matrix, in most cases, only one GCM member is downscaled. Leduc et al (2019); von Trentini et al (2019) use a SMILE ensemble to assess the role of the internal variability in future projection. Both studies confirm that the internal variability is more substantial at local scales, and CORDEX ensembles should better explore this aspect. The constitution of high resolution (2 to 25 km) large ensembles accounting for all sources of uncertainties is one of the biggest challenges of the regional climate community (Giorgi, 2019; Rummukainen et al, 2015).

Various studies proposed statistical approaches to fill up Scenario-GCM-RCM matrices (Déqué et al, 2007, 2012; Christensen and Kjellström, 2022; Evin et al, 2019). These methods are based on variance decomposition (ANOVA) and allow us to estimate the proportion of the different sources of uncertainties. Among the others, Evin et al (2019, 2021) propose a Bayesian ANOVA model to deal with an unbalanced matrix. Figure 1.11 illustrates the variance decomposition results for temperature and precipitation in winter and summer over three large European regions. Compared with the fractional uncertainty decomposition from Figure 1.8, the RCM is added, as well as the interaction terms between GCM, RCM and RCP scenarios. The results are comparable with the GCM uncertainty analysis, the internal variability role is decreasing through time, conversely to the scenario uncertainty. Regarding specific variables/seasons/regions, the uncertainty brought by the RCM can be equivalent to the one from the GCM. The internal variability seems to be more critical when looking at smaller regions. It is confirmed in Evin et al (2019), which conducts the same analysis at very local scales, in 4 French mountain massifs of the Pyrenees and Alps. The results show that internal variability is the largest source of uncertainty for temperature until the second half of the century, while for
Orography of the EURO-CORDEX domain at typical GCM, RCM and
CPM resolution

<figure><img src="image_14.png" /><figcaption>**Figure 1.12:** Illustration of typical resolutions for GCM, RCM and the very high resolution CPMs.</figcaption></figure>

precipitation, it vastly overcomes all other sources all along the period.

A well-known example of RCM added value: Extreme precipitation over
the Alps

The better representation of complex terrains such as reliefs, coastlines or is-
lands due to the higher resolution implies a better simulation of the local climate as
these are important drivers, and their interaction with the atmospheric circulation
is better represented. It is the first source of added value that benefits all variables.
Careto et al (2022) showed, for example, the added value of RCM over the Iberian
Peninsula regarding temperature. Nevertheless, the substantial and most recognised
added value of RCM concerns the representation spatial and temporal behaviour of
precipitation.

Multiple studies have shown that RCMs with higher resolution represent better the extreme precipitations over the Alpine region. Torma et al (2015) and Fantini et al (2018) notably showed that the CORDEX simulations containing the Alps recreate well the entire distribution of precipitation, both studies based on different scores of distribution similarity (Kolmogorov-Smirnoff and Kullback-Leiber). Figure 1.13 (a) illustrates the better representation of the mean daily rainfall by RCM. It is not shown here, but Torma et al (2015) shows that the added value of RCM is still visible at the GCM scale when interpolating (upscaling) the high-resolution simulation on the GCM grid. The second part of the figure shows that this better representation can modify the local climate change response signal. Indeed Giorgi et al (2016) showed that the RCM ensemble agreed on an increase of precipitation over the top of the Alps while the GCM ensemble simulated a strong uniform decrease of prediction. The authors found the RCMs projection more accurate, as they match with recent trends, and the underlying process checked with other variables is credible. The evaluation of climate models and RCMs in a climate change context cannot be empirical. Here is the strength of RCMs, the complete analysis of underlying physical processes is possible, and the credibility of the scenario can
<figure><img src="image_15.png" /><figcaption>**Figure 1.13:** (a) Summer mean precipitation over the Alps, comparison of RCM and GCM ensembles with versus the observed mean, from Torma et al (2015); Giorgi (2019).
(b) Projection of summer precipitation changes in GCM and RCM ensembles, extracted from Giorgi et al (2016); Giorgi (2019).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 1.14:** Extreme precipitation distribution. The plots show the returned CDF zooming on the extreme quantiles. It compares a CPM (AROME, 2.5km) and its driving RCM (ALADIN63, 12km) with two observations set, at 8km resolution and at 1km.</figcaption></figure>

be addressed. Besides, like for GCMs, the agreement of multiple GCM-RCM pairs bring robustness to the result.

Recent literature (Fumière et al, 2020; Berthou et al, 2020; Ban et al, 2021) have shown that the use of higher resolution CPMs largely improves the representation of extreme precipitations compared with 12km resolution RCMs. Figure 1.14, extracted from Fumière et al (2020) compares the distribution of extreme precipitation for a CPM model and its driving RCM versus two sets of observation. We can see that the CPM reproduces until the highest quantiles of the high-resolution observations. Caillaud et al (2021) shows that CPMs can recreate realistic objects of extreme precipitation thanks to a good representation of hourly rainfall. The time evolution and movements, the intensity and the frequency of the objects are well represented by the evaluated model in the evaluation running mode. The future projections indicate an intensification of such extreme events, but such findings need to be confirmed by ensemble studies.

## More on the RCM vs GCM large scales inconsistencies

Laprise et al (2008) propose an original evaluation of regional climate models following a set of 4 fundamental tenets that represent the initial paradigm of dynamical downscaling:

* Tenet 1: RCMs are capable of generating small-scale features absent in the driving fields supplied as lateral boundary conditions;

* Tenet 2: The small scales that are generated have the appropriate amplitudes and climate statistics;
* Tenet 3: The generated small scales accurately represent those that would be present in the driving data if it were not limited by resolution;

* Tenet 4: In performing dynamical downscaling, RCM generated small scales are uniquely defined for a given set of LBC.

The study easily validates the two first tenets, while the two last are mostly invalidated. Indeed, as mentioned earlier, the RCM is not fully deterministic and tends to derive from the large driving scale creating its own climatology. Moreover, several studies (Sanchez-Gomez et al, 2009; Lucas-Picher et al, 2008) have shown the impact of the initial conditions on the internal variability of the model, confirming the non-verification of tenet 4. Nevertheless, in climate studies, the deterministic property is not necessary, as only the statistical representation of the local climate is really of interest. Laprise et al (2008) proposes a last tenet which is declined in 3 versions as there is little consensus in the community. It concerns the modification of large scales by the RCM :

* Tenet 5a: The large scales are unaffected within the RCM domain;

* Tenet 5b: The large scales may be improved owing to reduced truncation and explicit treatment of some mesoscale processes with increased resolution within the RCM domain;

* Tenet 5c: The scales larger than or comparable to the RCM domain are degraded because the limited domain is too small to handle these adequately.

If the study invalidates tenet 5a, it is more difficult to assess whether the RCM improves or degrades the large scales. This is still a hot topic in the regional climate modelling community. The debate is intrinsically linked with the spectral nudging debate mentioned earlier. Various recent studies have shown that the Euro-Cordex ensemble largely underestimates the summer European warming projected by the CMIP5 models (Boé et al, 2020; Ribes et al, 2022). Here again, the opinion diverges on the soundness of this result. Some studies argue that the RCMs projection is more reliable than GCM because they show less bias in the historical period or because of the better representation of local physical processes impacting the large scales (Bartók et al, 2017; Sørland et al, 2018). On the other hand, other studies point out the absence of time-varying anthropogenic aerosols (Boé et al, 2020; Gutiérrez et al, 2020) and CO2 effect on vegetation (Schwingshackl et al, 2019; Boé, 2021) in most EURO-CORDEX simulations which would go in the line of CMIP5 models. Lastly, with a dedicated modelling experiment, (Taranu et al, 2022) concludes that the aerosols play a significant role in these discrepancies but also highlight the inconsistencies in physical parametrisations in GCM-RCM pairs.

### 1.3.2 Statistical Downscaling

The alternative to dynamical downscaling and RCMs relies on statistical approaches. Empirical/Statistical Downscaling (ESD) aims to estimate the link between high resolution observations (predictands, target, Y) and large scale climate variables (predictors X) (Maraun and Widmann, 2018; Maraun et al, 2010):

$$ \mathbb{E}(Y|X) = f_{\Omega}(X) \qquad (1.1) $$
with $\mathbb{E}$ the common expectation and $\Omega$ represents a set of unknown parameters estimated during the training or calibration phase. These methods focus mainly on surface variables, such as temperature, wind or precipitation, essential for impact studies. The model $f$ will vary according to the variable to downscale. The set of parameters $\Omega$ represents the local characteristics influencing the climate and will depend highly on the location. This estimated relationship $f$ is then applied to low resolution global climate model outputs to reconstruct high-resolution climate information. Over the last three decades, many methods have been proposed. Maraun et al (2010) regroups ESD methods into three classes that we will detail below.

## Perfect Prognosis

The perfect prognosis (PP) statistical downscaling methods aim to link a set of observed large scale predictors ($X = X_1, X_2, ..$) and observed local target variable $Y$ at time $t$. Eq. 1.1 becomes then:

$$ \mathbb{E}(Y_t|X_t) = f_{\Omega}(X_t) + \epsilon_t \quad (1.2) $$

With $\epsilon$, the error term represents what is not explained by $X$. The development of PP approaches includes three steps: selecting and preparing inputs (1), obtaining a good<fnref n="9" /> dataset for $Y$ (3), and choosing the statistical model (2). Accordingly, it is necessary to calibrate the model, i.e. estimate the parameters $\Omega$ with the observed data. Each $X_t$ has a unique corresponding $Y_t$, so the calibration is done element-wise. The observed $Y_t$ is generally a meteorological station allowing long period records. Recently, high-resolution gridded observational products have also been used since their period extent starts to be long enough (Baño-Medina et al, 2020). The predictors $X_t$ generally come from reanalysis that ensures homogeneous and continuous data.

Predictors are chosen to describe the general atmospheric circulation. Temperature, altitude of given pressure levels (geopotential), humidity or wind components are the most commonly used predictors (Gutiérrez et al, 2019). They can be considered at the surface level or higher altitudes. Multiple studies have shown that the best set of predictors highly depends on the season and the location (Timbal et al, 2009; Lemus-Canovas and Brands, 2020; Erlandsen et al, 2020). Moreover, those predictors are high-dimensional 2D fields, and it is necessary to reduce their dimension to suit most statistical models. Because of the high correlation between neighbour points and between variables, the application of variance decomposition methods such as Principal Components Analysis (PCA) helps to reduce the dimension drastically while keeping a maximum of information. Other methods taking into account the predictand can be applied for the same purpose. The high spatial variability of climate when it comes to very local scales, especially for precipitation, implies building very local models or choosing common predictors for larger target domain, which can degrade performance.

A wide range of statistical models has been proposed to estimate that relationship. To keep it simple, we can regroup them through two big families: the transfer

<footnote n="9">⁹i.e. long, homogeneous</footnote>
functions (TF) and analogues approach (A) (Gutiérrez et al, 2019; Vaittinada Ayar et al, 2015). The transfer functions methods generally rely on regression functions. Multiple Linear function have been widely used to model temperature (Huth et al, 2015; Hertig et al, 2013; Gutiérrez et al, 2013) with various predictor selection techniques. The highly skewed nature of precipitation distribution requires the use of Generalized Linear Models (GLM) as they allow Y to follow a non-Gaussian distribution (San-Martín et al, 2017). The linear model estimates a transformation of the expected outcome (see Eq 1.3). More complex regressions such as General Additive Models (GAMs) that generalize GLMs to nonlinear models, Vector GLM (Vrac et al, 2007a) or non-parametric regressions are existing alternatives.

$$g(\mathbb{E}(Y|X)) = \omega_1 X_1 + \omega_2 X_2 + \dots \qquad (1.3)$$

with g the transformation function and $\omega_i \in \Omega$ the model weights.

Recent studies explored the application of machine learning algorithms, notably deep neural network architecture (Vandal et al, 2019; Baño-Medina et al, 2020, 2021; Höhlein et al, 2020). The big advantage of such methods is their ability to identify the appropriate predictors (Baño-Medina et al, 2020) and to deal with high-dimensional inputs and outputs. Indeed, most TF methods are univariate (single-site, single variable), while deep neural networks can easily produce maps and can explore more than one variable. Neural networks are an essential part of this work, and their functioning will be detailed in chapter 2, later in the manuscript. Most methods do not account for $\epsilon$, making them fully deterministic, but some are stochastic since they propose an estimation of the variance or distribution parameters (Cannon, 2008; Baño-Medina et al, 2020)

On the other hand, analog methods use the predictors $X_t$ to find similar (analog) atmospheric situations $A(X_t)$ in the calibration set, such that the $\mathbb{E}(Y_t|X_t)$ is a weighted combinations of $Y_{A(X_t)}$. Those methods have primarily been used in weather prediction (Lorenz, 1969). Different variations on the way to evaluate the distance between predictors or the way to combine the predictions exist (Dayon et al, 2015; Maurer et al, 2010; Gutiérrez et al, 2013; Vaittinada Ayar et al, 2015).

## Model Output Statistics (MOS)

MOS approaches aim to link model outputs to observed station records directly. The idea is to find the function that transforms the model output distribution into the station data distribution to fit the local climatology. MOS methods do not imply predictor selection. The target variable low resolution model output is the only predictor and it is taken from the closest grid point. Thus, MOS approaches are quite simple to use. However, the relation is valid only for the calibration model, as the biases are model dependent. Most of these methods are based on quantile mapping algorithms, and a wide range of variants have been developed (Déqué et al, 2007; Michelangeli et al, 2009; Štěpánek et al, 2016; Themeßl et al, 2012). For precipitation, some methods comport two steps: the prediction of occurrence first and then the amounts on wet days. This is, for example, the case of Volosciuk et al (2017), which uses logistic regression to estimate the occurrence or not of rain and a GLM to fit the parameters of a gamma distribution. The great majority of
<figure><img src="image_17.png" /><figcaption>**Figure 1.15:** Example of weather generator (from Maraun et al, 2010; Vrac and Naveau, 2007) based on non-homogeneous Markov chains. The transition from one state to another depends on the past state (hidden Markov model) and the atmospheric conditions. Moreover, the amount of rain on a wet day is generated from a distribution dependent on the atmospheric conditions.</figcaption></figure>

MOS approaches are single site and single variable. Nevertheless, MOS methods are similar in many points to Bias Correction techniques. Recent developments of multivariate BC techniques (Robin et al, 2019; Vrac, 2018; François et al, 2020, 2021) could eventually be applied in statistical downscaling framework.

## Weather Generator (WG)

Weather Generators are defined as stochastic models of meteorological variables that explicitly model their marginal distribution and temporal dependence (Maraun and Widmann, 2018). They regroup various models, including Markov chain or auto-regressive models. The model parameters are fitted thanks to observed time series. For instance, first-order Markov chains have been widely used to simulate rain occurrence (Richardson, 1981; Katz, 1977; Wilby et al, 1998; Dubrovský et al, 2004) and are coupled with gamma distribution to model the precipitation amount. Alternatively to the gamma distribution, some studies have used a mixture of exponential distribution (Keller et al, 2017).

In the downscaling framework, statistical changes from one period to the others are evaluated in a climate model at the grid point corresponding to the observed series. Those changes are then applied to the statistical model (Kilsby et al, 2007). Another way is to include covariates from the climate model at different steps of the statistical model. Figure 1.15 (from Maraun et al, 2010; Vrac and Naveau, 2007) illustrates such weather generator which can also be considered as a perfect prognosis method.
**Evaluation and intercomparison studies**

The fifth assessment report of the IPCC concluded that it is difficult to provide
a general assessment of ESD approaches performances because, among other things,
of the wide variety of approaches. In most studies mentioned above, ESD methods
are evaluated independently on their ability to reproduce a hidden set of observed
series. Some intercomparison studies exist, but they mainly focus on a particular
family of methods, such as MOS approaches in Gudmundsson et al (2012); Gut-
mann et al (2014). For this purpose, the VALUE framework (Maraun et al, 2015),
proposes a large intercomparison exercise over Europe, where different teams can
propose their downscaling methods. The VALUE project includes three steps. A
first perfect predictor experiment evaluates the methods versus observations using
reanalysis predictors. Then it plans the inter-comparison of the methods with GCM
historical (Experiment II) and future (Experiment III) simulations. At the date
where this manuscript is written, only the conclusions of the first experiment have
been published (Gutiérrez et al, 2019; Hertig et al, 2019; Soares et al, 2019; Maraun
et al, 2019a,b). Recently, Vaittinada Ayar et al (2015) proposed also an intercom-
parison of downscaling methods including both RCMs and ESD over Europe and
Bettolli et al (2021) one over South America in the CORDEX framework.

The general result of ESD approaches in perfect predictor experiments is that
a well-chosen and constructed model can fit the original data correctly. MOS ap-
proaches are generally simple to set up and provide good results. The performance
of PP methods relies a lot on predictor choice, and a very local model is preferable
to a large one. WG can recreate the correct statistics properties of the original
series but are unable (by construction) to create appropriate temporal and spatial
variability with, for example, the case of dry or hot/cold spells. More generally, the
representation of extreme values is more difficult (Hertig et al, 2019), even if some
methods perform better than others (Vaittinada Ayar et al, 2015). The critical point
is that no method seems to outperform another in all aspects. The second point
is that ESD methods can represent realistic high-resolution surface variables in the
present climate. Compared to RCMs, they often have fewer biases with respect
to observations, which seems logical since they use them for calibration. Vandal
et al (2019) or Baño-Medina et al (2020) propose an intercomparison of methods,
including recent machine learning algorithms. Neural Networks showed outstanding
accuracy in (1) dealing with gridded inputs and outputs and (2) reproducing well
the statistical properties of the target variables.

Because of the evident lack of future observations, it is difficult to assess whether
the estimated relationship between large and local scales holds in a future or different
climate. Indeed, applying ESD approaches to climate model projections relies on the
strong assumption that the relationship is stationary. In order to assess whether this
hypothesis holds or not, some studies have used an RCM as observations, following
a perfect model framework (Vrac et al, 2007b; Gutiérrez et al, 2013; Dayon et al,
2015; Lanzante et al, 2018; Erlandsen et al, 2020). *f* is estimated in the historic run,
and tested on future projections. All those studies conclude that ESD methods are
performing better in the historical period (where the relationship is learnt) than in
future, especially regarding extremes (Vrac et al, 2007b). Lanzante et al (2018) found
that the ESD method applied could not reproduce the RCM changes over mountains
area and along the littorals, regions where the RCM bring added value thanks to its high resolution. Gutiérrez et al (2013) also concludes that surface predictors are preferable because of a changing relation between altitude and surface temperature. Finally, Erlandsen et al (2020) finds that the calibration period and length are key factors for the method’s accuracy. Besides the stationarity assumption, which seems not to hold, this can be another substantial limitation for ESDs methods. If there are numerous long-term records of observations in Europe and North America, it is not the case for most places worldwide (Bettolli et al, 2021).

1.4 Objectives and scientific questions of this study.

We have seen in this introductory chapter that there is a strong need for reliable
information about local climate change to conduct impact studies and prepare adaptive strategies and policies. “Reliable” implies a proper exploration of the different sources of uncertainties to assess what is very likely to happen or what is among the possible. It is also essential to explore the extreme worst-case scenarios among the possibles, as they imply huge societal costs.

Today the primary tool to study climate change is the family of global climate models and earth system models. However, their coarse resolution does not allow evaluation of the impacts at small scales. To answer this issue, two downscaling communities have emerged to fill the gap between local climate and GCM large scales. Empirical Statistical Downscaling methods present the decisive advantage of being relatively not expensive and quite easy to set up. However, they suffer some lack of trustability to reproduce future changes as they rely on the *stationarity assumption* (Dixon et al, 2016). Besides, the calibration requires long-term records of observations which are difficult to find in most places worldwide. On the other hand, dynamical downscaling relies on regional climate models that are recognised for their precise representation of local physical processes and simulate a realistic local climate. For instance, they have good abilities to simulate extreme events, dry/hot spells or heavy precipitations. However, their high expensive costs strongly limit the exploration of the uncertainties.

In this context, this work aims to answer the strong challenge of producing re-
liable information on local climate change by reconciling these two fields. Indeed,
the strengths and limitations of each family of methods match the others pretty
well. This work aims to expand the number of high-resolution RCM simulations
regarding some specific variables. The idea is to learn, with a machine learning
algorithm, the downscaling function included in the RCM to apply it to other GCM
simulations. According to the definition given in Reichstein et al (2019)<fnref n="10" />, we aim
to build an Emulator for the RCM downscaling function. We will name this novel
approach RCM-Emulator. Another way of seeing this work’s objective is to build
an ESD approach that uses RCM simulations to be calibrated. The advantage of
this approach is the possibility of learning the downscaling relationship in different

<footnote n="10">¹⁰Emulation: Emulation of the full (or specific parts of) a physical model can be useful for computational efficiency and tractability reasons. Machine learning emulators once trained can achieve orders of magnitude faster simulations than the original physical model without sacrificing significant accuracy</footnote>
climates (including warmer climates) and regions with no observation records. The RCM-emulator proposed here could be categorised as a Perfect Prognosis approach as it aims to learn the relationship between low-resolution large-scale predictors and high resolution surface variables.

Other studies have already combined ESD techniques with RCMs, notably in 2-step downscaling approaches where predictors are taken from RCMs simulations (Vrac et al, 2012; Pryor and Barthelmie, 2014). This is particularly the case for MOS approaches initially designed for this purpose. Chadwick et al (2011) is probably the first study to propose combining ESD and RCM similarly to this work. The authors use neural networks to link GCM and RCM outputs. The study showed promising results and notably concluded the importance of learning in future climates. Walton et al (2015) use RCM simulations to estimate the monthly pattern of climate change over California and then create a small ensemble by emulating five other GCMs. However, it used a complex and domain-specific method to select the predictors, and the model was calibrated concisely (10 years in the past and future). The RCM-Emulator proposed in this work relies on deep neural network algorithm as they have already shown promising results in a perfect prognosis context. Recently, multiple studies have proposed similar approaches also based on neural network architectures (Wang et al, 2021; Serifi et al, 2021; Babaousmail et al, 2021). The methods proposed are closer to a MOS approach as they aim to recreate the high resolution temperature or precipitation from the low resolution of the same variables.

Hereafter we introduce the scientific questions that will help us to define, develop and validate the novel hybrid downscaling tool we propose here. A complete understanding of the strengths and weaknesses of the tool is of primary importance in assessing its reliability.

1. **Is the RCM-Emulator feasible?**

The first essential step is to define the concept of the emulator properly. It implies defining the exact objective and the adapted framework to isolate and estimate the downscaling function carried out in a Regional Climate Model. Moreover, the choice of the statistical method used to estimate it is also part of the conceptualisation of the emulator.

2. **How reliable is the emulator?**

The second important aspect in developing a novel approach is guaranteeing its trustability. For both downscaling families, numerous studies have been dedicated to proving the reliability of the downscaling method. It implies finding the adapted frameworks to study the downscaling tool's abilities in different steps. It is notably necessary to evaluate the emulator production to reproduce the known added value of RCM. Moreover, the proper identification of strengths and weaknesses is essential to give confidence in the tool.

3. **How efficient?**
Efficiency is a crucial attribute of the emulator. The tool aims to tackle the high costs of RCMs and expand the number of downscaled simulations. The efficiency objective must be kept in mind throughout the emulator’s development and validation.

4. How transferable?

Finally, the transferability of the emulator to different levels needs to be as-
sessed. The emulator will be defined in a certain framework, and evaluating
how it reacts when moving out of this structure is important. The transfer-
ability can be considered to other sources of input, to other variables, to other
climates or different downscaling problems. Assessing the transferability of the
emulator is necessary to define how and to what problems it can be applied.

The RCM we will emulate in this study is ALADIN63 (Nabat et al, 2020). It is the model developed at CNRM, the research centre associated with the French weather service and CNRS. The ALADIN63 physical scheme is the same than the CNRM-CM family of models also developed at CNRM (in collaboration with CER-FACS). ALADIN63 has been firstly used to downscale CNRM-CM5 (Voldoire et al, 2013), and then three others GCM: MPI-ESM-LR, developed at the Max Plank Institute (Giorgetta et al, 2013), NorESM1-M, developed at the Norvegian Climate Center (Bentsen et al, 2013) and the HadGEM2-ES, the Met Office climate model (Collins et al, 2011). This small ensemble is a good opportunity to build and evaluate the RCM-Emulator proposed in this study. The better consistency between ALADIN63 and CNRM-CM5 offers the best case for the first tests.

This manuscript is organised into five chapters after this introduction. Chapter 2 is entirely dedicated to the machine learning algorithm we chose for the emulator. As mentioned earlier, neural networks have shown accurate results in statistical downscaling. We will explore their functioning and analyse why those methods are well adapted for climate data and especially downscaling problems. Chapter 3 presents a feasibility study of the emulator. It aims to describe the concept and set the framework to train and evaluate the RCM-emulator for temperature. Chapters 4 and 5 explore then two transferability aspects of the emulator. Chapter 4 assesses its robustness to other input sources, while Chapter 5 explores the application to precipitation which is a more challenging variable. Finally, the last chapter derives some conclusions regarding this work’s achievements and introduces some perspectives that emerge regarding the development and use of RCM-Emulators.
# Chapter 2

## Why deep neural networks are a credible strategy to build the RCM-emulator ?

This chapter aims to understand the general functioning of artificial neural networks and why they are well adapted to build the RCM-Emulator. We will review the critical aspects of neural networks which led to the recent boom of deep learning approaches. I will specifically focus on convolutional neural networks and the specific UNet architecture that I chose for the Emulator because of its ability to deal with images. Finally, I will review some examples of deep learning applications in climate and weather sciences.

### 2.1 A neural network

Introduced by McCulloch and Pitts (1943), the artificial neuron is a mathematical representation of the biological neuron. Eq. 2.1 describes the behaviour of such a neuron. It receives an input $X$ of dimension $n$ and returns the evaluation $y$ of dimension 1 through the activation function $g$ of the weighted combination of the inputs. [$\omega_1, \dots, \omega_n$] is the weights vector associated with the inputs. $\omega_0$ is a constant term called the neuron bias<fnref n="1" />. The activation function $g$ is a non-linear function determining the ‘type’ of neuron. The different activation functions commonly used are shown in Table 2.1.

$$ 
\begin{gathered}
y = f_{\Omega}(x_1, \dots, x_n) = g(\omega_0 + \sum_{i=1}^{n} \omega_i x_i) \qquad (2.1) \\
\text{with } \Omega = \omega_0, \dots, \omega_n
\end{gathered}
$$

The simple neuron model can be seen as a simple linear binary classifier. It is the perceptron model, introduced by Rosenblatt (1958). The author proposed an algorithm to find the suitable weight vector that correctly classifies each input $X$. The weights’ vector defines the hyperplane’s coordinates that divide the input space into two distinct regions. The algorithm evaluates a learning set of input-outputs pairs

<footnote n="1">¹An alternative expression consists to set $x_0 = 1$ and then $y = g(\sum_{i=0}^{n} \omega_i x_i)$.</footnote>
<figure><img src="image_18.png" /><figcaption>Figure 2.1: Neural Network diagram</figcaption></figure>

$(X_t, Y_t)_{1 \le t \le N}$ and updates at each iteration the weight vector in case of misclassification. The algorithm's convergence has been proven in Novikoff (1962). Rosenblatt already imagined that two connected neurons could solve non-linear classification problems, but it was computationally too expensive to be implemented.

A neural network is then a set of neurons organised through layers. A network includes three types of layers illustrated in Figure 2.1: the input and output layers and some hidden layers between them. Hidden layers perform non-linear transformations of the inputs to match the targeted output. A hidden layer includes multiple neurons that are not inter-connected. The organisation of the neurons inside a layer depends on the type of layers and the architecture chosen. For example, in a fully connected network, all neurons from a given layer are connected to all neurons from the previous one. Other architectures exist such as Convolutional Neural Network (CNN) that I will present later on. The idea behind a neural network is to mimic the brain; each neuron can only perform a simple action, but the whole network can solve highly complex non-linear problems.

A neural network is then a set of parametric functions $F_\theta$ which links some inputs $X_t \in X$ and outputs $Y_t \in Y$. $\Theta$ represents the parameter space that regroups all possible values of $\theta$, i.e. all the weights $\omega_{i,j}^k$ for the $i$-th input of the $j$-th neuron of layer $k$ plus the bias $\omega_{0,j}^k$. An error (or loss) function is associated with the neural network. Given a training set $\mathcal{D} = \{(X_t, Y_t), t \in [1, N]\}$, the error function $L(D, \theta)$ evaluates the distance between the target outputs $Y_t$ and the actual network outputs $\hat{Y}_t$ for the input $X_t$. The mean square error (Eq. 2.2) is the most common loss function for regression problems. Training the network involves finding the best set of parameters $\theta^*$ to minimise the training set's error:

$$L_{MSE}(D, \theta) = \frac{1}{N} \sum_{i \in \mathcal{D}} (Y_i - \hat{Y}_i)^2 \qquad (2.2)$$
<table><caption>**Table 2.1:** Commonly used activation functions</caption><thead><tr><th>NAME</th><th>EXPRESSION</th></tr></thead><tbody><tr><td>Identity</td><td>g(x) = x</td></tr><tr><td>Heaviside</td><td>g<sub>a</sub>(x) = <span style="font-size: 2em; vertical-align: middle;">{</span> 0 if x &lt; a<br>1 otherwise <span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>sign</td><td>g<sub>a</sub>(x) = <span style="font-size: 2em; vertical-align: middle;">{</span> -1 if x &lt; a<br>1 otherwise <span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>sigmoid</td><td>g(x) = <div style="display: inline-block; vertical-align: middle; text-align: center;"><span style="display: block; border-bottom: 1px solid black; padding: 0 0.2em;">1</span><span style="display: block; padding: 0 0.2em;">1 + exp(-x)</span></div></td></tr><tr><td>tanh</td><td>g(x) = <div style="display: inline-block; vertical-align: middle; text-align: center;"><span style="display: block; border-bottom: 1px solid black; padding: 0 0.2em;">exp(x) - exp(-x)</span><span style="display: block; padding: 0 0.2em;">exp(x) + exp(-x)</span></div></td></tr><tr><td>Rectified Linear Unit (ReLU)</td><td>g(x) = <span style="font-size: 2em; vertical-align: middle;">{</span> 0 if x &lt; 0<br>x otherwise <span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>Leaky RELU</td><td>g<sub>α</sub>(x) = <span style="font-size: 2em; vertical-align: middle;">{</span> αx if x &lt; 0<br>x otherwise <span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>softplus</td><td>g(x) = ln(1 + exp(x))</td></tr></tbody></table>

Therefore, the learned network corresponds to:

$$F_{\theta^*} = \arg \min_{F_{\theta, \theta \in \Theta}} L(\mathcal{D}, \theta) \qquad (2.3)$$

To simplify, we include the bias $\omega_0$ in the weight vector $\omega$.

## 2.2 Training the network: the backpropagation algorithm

The next significant improvement in neural networks is the backpropagation algorithm introduced by Rumelhart and Hinton (1986). The big challenge when calibrating a neural network is to find the proper weights and biases for hidden layers while each neuron does not have a defined target outcome. Each neuron must be calibrated following the entire network error, implying that each neuron error function depends on the parameters from all previous and following layers. A direct update of all parameters, like in the perceptron model, is unimaginable because of the computations it implies. The backpropagation algorithm is an efficient optimisation scheme based on the gradient descent algorithm to address this issue. I briefly describe the gradient descent below.
<table><caption>**Table 2.2:** Notations for the backpropagation algorithm</caption><thead><tr><th colspan="2">Notations</th></tr></thead><tbody><tr><td>X &isin; &#8477;<sup>n</sup></td><td>Set of inputs to the network</td></tr><tr><td>Y &isin; &#8477;</td><td>1-D (for this example) target variable</td></tr><tr><td>&Ycirc;<sub>t</sub></td><td>Output of the network for sample t</td></tr><tr><td>D = {(X<sub>t</sub>, Y<sub>t</sub>), t &isin; [1, N]}</td><td>Training set</td></tr><tr><td>L</td><td>Loss function</td></tr><tr><td>&theta;</td><td>Set of parameter for the entire network</td></tr><tr><td>&Omega;</td><td>Set of weights for one neuron</td></tr><tr><td>k</td><td>Indicator for the layer, the last layer is called z</td></tr><tr><td>j</td><td>Indicator for the neuron in layer k</td></tr><tr><td>i</td><td>Indicator for the input of the jth neuron in layer k</td></tr><tr><td>&omega;<sub>i</sub>, j<sup>k</sup></td><td>weight for the input i of the jth neuron in layer k</td></tr><tr><td>o<sup>k</sup><sub>j</sub></td><td>Output of the jth neuron in layer k</td></tr><tr><td>O<sup>k</sup></td><td>Number of outputs for layer k</td></tr><tr><td>a<sup>k</sup><sub>j</sub></td><td>Activation of the jth neuron in layer k</td></tr><tr><td>g, g<sub>z</sub></td><td>Activation function, (z for the activation function of the last layer)</td></tr><tr><td>&ell;</td><td>loss function output for one given sample (X*, Y*)</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Figure 2.2:** Illustration of a neuron with the notation for Section 2.2</figcaption></figure>
The gradient descent is an iterative algorithm that finds a local minimum of a differentiable convex function. The gradient at position $p$ of a $n$-dimensional function $h$, written $\nabla h(p)$ is the vector of its $n$ partial derivative evaluated at this point. It gives the direction of the “fastest” increase of $h$ at a given point $p$. In order to reach a local minimum, the algorithm evaluates the function gradient at a starting point and moves in the gradient’s opposite direction following a pre-defined step-size called “learning rate”. By iterating this process, one should be able to find a local minimum of the given function following Eq 2.4.

$$p_{t+1} = p_t - \lambda \nabla h(p_t) \qquad (2.4)$$

In the case of neural networks, the aim is to minimise the loss function $L$ according to the set of parameters $\theta$. The backpropagation algorithm proposes an efficient way to update all the weights and biases of the network by simplifying the gradient descent computation. The idea is to evaluate the gradient in each weight $w_{i,j}^k$. Here we consider the loss function to be MSE which can be decomposed into a sum of individual errors (Eq. 2.5). It allows us to simplify the explanation here and to look at the computation for a single couple $(X^*, Y^*)$ and we write the loss function evaluated at this point $\ell$ (Eq. 2.6). The notations for this part are summarized in Table 2.2 and illustrated for one neuron in Figure 2.2.

$$\frac{\partial L(\mathcal{D}, \theta)}{\partial \omega_{i,j}^k} = \frac{1}{N} \sum_{t \in \mathcal{D}} \frac{\partial L((X_t, Y_t), \theta)}{\partial \omega_{i,j}^k} \qquad (2.5)$$

$$\ell = L((X^*, Y^*), \theta) \qquad (2.6)$$

First of all we can re-write the partial derivative using the chain rule and defining the activation $a_j^k$ that is the weighted sum of the neuron inputs plus the bias.

$$\frac{\partial \ell}{\partial \omega_{i,j}^k} = \frac{\partial \ell}{\partial a_j^k} \frac{\partial a_j^k}{\partial \omega_{i,j}^k} \qquad (2.7)$$

We set the first term as $\delta_j^k = \frac{\partial \ell}{\partial a_j^k}$. It somehow represents the contribution of the $j$th neuron of layer $k$ to the total error $(\hat{Y}^* - Y^*)$. The second term $\frac{\partial a_j^k}{\partial \omega_{i,j}^k}$ is the derivative of a sum. It is simply equal to the $i$th input of the neuron, which is the output $o$ from the $i$th neuron of the previous layer<fnref n="2" /> (Eq. 2.8). We define $O^k$ as the total number of outputs of layer $k$.

$$\frac{\partial a_{i,j}^k}{\partial \omega_{i,j}^k} = \frac{\partial \sum_{r=0}^{O^k-1} \omega_{r,j}^k o_r^{k-1}}{\partial \omega_{i,j}^k} = o_i^{k-1} \qquad (2.8)$$

The missing terms to evaluate are the $\delta_j^k$. For the last layer $z$, we consider having a single output model so $j = 1$, we have:

$$\delta_1^z = \frac{\partial \ell}{\partial a_1^z} = \frac{\partial \frac{1}{2}(g_z(a_1^z) - Y^*)^2}{\partial a_1^z} = (g_z(a_1^z) - Y^*)g'_z(a_1^z) = (\hat{Y}^* - Y^*)g'_z(a_1^z) \qquad (2.9)$$

<footnote n="2">²We assume here to be in a fully connected network where all neurons from a layer are connected to all neurons from the previous and subsequent layers.</footnote>
since $\hat{Y}^* = g_z(a_1^z)$ with $g_z$ the activation function of the last neuron.

For the hidden layer, we can write $\delta_j^k$ according to $\delta_j^{k+1}$ thanks to the backpropagation formula :

$$ 
\begin{aligned}
\delta_j^k &= \frac{\partial \ell}{\partial a_j^k} = \sum_{r=1}^{O^{k+1}} \frac{\partial \ell}{\partial a_r^{k+1}} \frac{\partial a_r^{k+1}}{\partial a_j^k} \\
&= \sum_{r=1}^{O^{k+1}} \delta_r^{k+1} \frac{\partial a_r^{k+1}}{\partial a_j^k} \\
&= \sum_{r=1}^{O^{k+1}} \delta_r^{k+1} \left( \frac{\partial}{\partial a_j^k} \sum_{s=1}^{O^k} \omega_{s,r}^{k+1} g(a_s^k) \right)
\end{aligned}
\qquad (2.10) $$

$$ \text{Finally, } \delta_j^k = g'(a_j^k) \sum_{r=1}^{O^{k+1}} w_{j,r}^{k+1} \delta_r^{k+1} $$

where $O^{k+1}$ is the number of neurons from the next layer<fnref n="3" /> and $w_{j,r}^{k+1}$ are the weights given to the output of the j-th neuron from layer k in layer k + 1. The formula gives that in each layer, the partial derivatives for each neuron's weights and biases can be evaluated from the following layers. Finally, the update formula for each weight, taking into account the entire training set $\mathcal{D}$ and a learning rate $\lambda$ is given by :

$$ \Delta \omega_{i,j}^k = -\lambda \frac{\partial L(\mathcal{D}, \theta)}{\partial \omega_{i,j}^k} = \frac{1}{N} \sum_{t=0}^{N} \frac{\partial \ell_t}{\partial \omega_{i,j}^k} = \frac{1}{N} \sum_{t=0}^{N} \delta_j^{t,k} o_{i,j}^{t,k-1} \qquad (2.11) $$

$$ \text{With } \delta_j^{t,k} = \begin{cases} (\hat{y} - y)g_z'(a_1^z) & \text{if } k \text{ is the output layer} \\ g'(a_j^k) \sum_{r=1}^{O^{k+1}} w_{j,r}^{k+1} \delta_r^{k+1} & \text{else} \end{cases} $$

The optimisation algorithm comports then 2 steps computed for each couple $(x_t, y_t)$. In the forward step, $x_t$ is given to the network and propagated through all layers to get the prediction. All weights and neurons' outputs are saved during this step. In the backward step, the $\delta$s are computed from the last layer and backpropagated. The partial derivative from layer k are computed after the ones from layer k + 1 following equation. After doing so for the entire training set, it is possible to calculate all partial derivatives and update $\theta$ 2.11.

The backpropagation algorithm was a huge improvement for neural networks since it was finally possible to train a neural network with multiple layers. From this point, there have been many developments with different layers or architectures, as we will see some examples in the next sections. The optimisation algorithm has also evolved to suit novel kinds of layers or to improve the convergence speed, but the principle of backpropagation is still the core of all variants. Most evolutions concern the gradient descent strategy. For example, the stochastic gradient descent works

<footnote n="3">³I use $r$ here to go through all input for a given layer, it ranges from 1 to leave room for $\omega_0$ which is the bias.</footnote>
on the same principle but updates $\theta$ using only one example from the training set and not the entire dataset at each iteration. In this work, I use the ADAM (Adaptive Moment Estimation) optimiser introduced in Kingma and Lei Ba (2015) and recognised for its efficiency and good performance. It is an adapted stochastic gradient descent based on estimating the gradient's first and second-order momentum to adapt the learning rate for each parameter.

## 2.3 The last milestones: the ReLu function and the computational resources

In neural networks development, we are missing two important points to fill the 25 years gap between the backpropagation algorithm in the mid-80s and the boom of neural networks in the early 2010s. The first point is related to the vanishing gradient. The sigmoid or tanh functions were the most used activation function until the study of Glorot et al (2011). They showed two strong advantages of using the ReLu function (see Table 2.1) for hidden layers. It is almost linear, making computation easier (so quicker) and avoiding any gradient vanishing effect. Indeed sigmoid or tanh, as they are functions from $\mathbf{R}$ to $[0, 1]$, tend to reduce the gradient after several steps in the backpropagation algorithm. The second advantage concerns the sparsity of the model, thanks to the negative part of the function. Indeed, ReLu can set to 0 some neurons bringing sparsity to the network, which Glorot et al (2011) and other studies found benefit for training the network.

The second point implies computational techniques and the new size of the available training dataset. In 2009, two major publications marked the beginning of neural network expansion. First, Deng et al (2009) created ImageNet, the current largest labeled image dataset. It constituted the largest training set, allowing to train deep neural networks (with numerous hidden layers). Lastly, Raina et al (2009) started to use Graphic Processor Unit (GPU) to train neural networks, which sped up the training of networks at an unprecedented rate (at least 10 times for the smallest networks and up to 70 times for the most time consuming ones).

## 2.4 Convolutional neural network (CNN)

The neural network presented until now and illustrated in Figure 2.1 is called a Fully Connected Network. It is composed only of Dense layers where each neuron is connected to all outputs from the previous layer. They have shown limitations in dealing with high-dimensional inputs, notably with images. Indeed as all neurons are connected, the number of parameters increases exponentially with the input size. Moreover, too many connections increase the risk of overfitting<fnref n="4" />. LeCun et al (1998), introduces LeNet, the first convolutional neural network specially designed to deal with images (visual recognition, classification, etc.). The originality of the network is the *convolutional layer* that I explain below.

<footnote n="4">⁴The neural network only learns “by heart” to reproduce the training set instead of estimating the function.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 2.3:** Convolution and Pooling principles, Goodfellow et al (2016).</figcaption></figure>

As indicated by its name, the convolutional layer performs a convolution be-
tween the input matrix (an image where each matrix coordinate is a pixel) and a
kernel matrix. The idea is to look at the image blockwise and extract information
from each block. The convolution process is illustrated in Figure 2.3. The kernel
moves through the entire image such that the convolution output is another matrix
of smaller size. This output matrix, called *Activation Map* on the figure, is then
passed to the activation function. The weights $\omega$ are then the values of the kernel
matrix. Instead of 16 weights to fit on the illustration of Figure 2.3 with a dense
layer, there are only 4 in this convolution layer. The number of parameters is defined
by the kernel size, not input size. In order to get various decompositions, multiple
kernels, called filters, can be passed through the input to isolate different aspects of
the input. Moreover, a colourful image is a 3D matrix, where the height and width
correspond to the image size, and the third dimension contains the three colour
channels (red, green, and blue). Various kernels are then passed on to the different
channels.

A strong advantage of a convolutional layer is that it uses the same patch of
parameters in the whole image with the reasonable assumption that if a given kernel
has been useful in a first position on the image, it is likely to be also useful somewhere
else. There are three main parameters to define in a convolutional layer :

* The depth: the number of filters to pass through the image. The size of the kernel also has to be defined.

* The stride: the number of pixels to move the kernel. On Figure 2.3 the stride is 1, with a stride equal to 2, the output would be a 2 × 2 matrix.

* The Zero-padding: it helps control the output size by adding 0 around the input borders.
The *Pooling layer* is the other specificity of CNNs. It is here to cut the spatial size of its input. It divides the input matrix into separated blocks of the same size and applies a simple function, such as the maximum or the mean, to each of them (cf. Figure 2.3). After a convolutional layer, a pooling window of size 2×2 downsizes each feature matrix by 2. There are two main advantages of pooling. First, it reduces the number of computations, making the network more efficient and avoiding overfitting. Second, it summarises the information extracted by the convolutional layer in a given region. It makes the model more robust to little variations.

Convolutional neural networks are a succession of convolutional and pooling layers for the extraction of features followed by a fully connected network which performs the classification or regression task as illustrated in Figure 2.4. Various CNNs have been developed since LeNet (LeCun et al, 1998): AlexNet introduced by Krizhevsky et al (2012) marked the beginning of neural network and CNNs lead for image recognition task or VGGnet introduced in Simonyan and Zisserman (2015) which is still widely used nowadays.

If CNNs suit image classification or object detection well, they do not suit the semantic segmentation problem. The aim is to assign each pixel to the different objects in the image. The target of such a problem is a mask with the input image size. Long et al (2015) propose a fully convolutional network to address this exercise. It takes the same base as a classical CNN but replaces the dense part with an upsampling layer that enlarges its input's spatial size. There exists multiple upsampling functions, one being the transpose convolution.

A very efficient and famous architecture for image segmentation is the UNet (Ronneberger et al, 2015), shown in Figure 2.5. The idea of this network is to divide the network task into two equal parts. In a first contractive path, the neural network looks like a classical CNN. It is composed of downsampling blocks, including convolutional layers and max pooling. At each step, the spatial dimensions are divided by two while the number of features doubles. It extracts the required features of the input image. The originality of the UNet resides in the expanding or upsampling path, which is almost symmetrical to the downsampling part (making the U-shape). It is composed of transposed convolutions to increase the spatial dimension and recreate the output size. Classical convolutions with multiple features allow the network to propagate information to the high-resolution layers. The downsampling blocks' outputs are concatenated with each upsampling layer, which should help the network to relocate each pixel correctly. The UNet architecture can be seen as an encoder/decoder network. Such networks aim to encode high dimensional and complex information in a latent space of reduced dimension while being able to reconstruct the original data with the decoder. In UNet, the aim is not to reduce dimension, but the important information is extracted and encoded in the space of a smaller dimension and then injected into the expanding part. The UNet will be an essential part of my work, notably because of its superior ability to extract information from a set of inputs and re-inject it to reconstruct a new image.
<figure><img src="image_21.png" /><figcaption>**Figure 2.4:** Illustration of CNN architecture, from Hoeser and Kuenzer (2020). (b) presents the entire network while (a) and (c) zoom on the beginning and end parts.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 2.5: Illustration of UNET architecture, from Ronneberger et al (2015).</figcaption></figure>

## 2.5 Other great families of architecture

On top of the family of convolutional neural networks, other families dedicated to specific tasks exist. I will briefly introduce two other noteworthy kinds of networks. The point here is not to detail their functioning but to understand the main idea behind these two families and for what they are designed.

First, Recurrent Neural Networks (RNNs) help to handle sequential data as time series or sentences. The idea behind RNN is to predict the output at time *t* thanks to the previous prediction and the corresponding input. Input and output are then considered sequentially. The most potent version of RNN is called LSTM and has shown good abilities to deal with long series. RNN can learn different types of relations as they can map one input to one output, one input to many outputs, many inputs to one output or many to many. Thus, they are widely used for text generation, time series forecasting or video anticipation.

The Generative Adversarial Network (GAN) is the other important family of neural networks. They consist in two neural networks that the user puts in competition. In general, they are used in the context of image generation. The first network is a generator which produces an image. The second network, called Discriminator, is trained to recognise if an image is true or fake, i.e. produced by the Generator. The two networks are then competing: the Generator tries to produce an image as close as possible to the “real” ones (i.e. such that the Discriminator classifies them as true ones). Simultaneously, the Discriminator keeps learning to differentiate between true and fake images. The network is trained when both Generator and Discriminator reach an equilibrium. Thus, the output of a GAN is one possible output among others. Neural networks generally produce what minimises the loss function, which can be a compromise of the less lousy output. A way to
<figure><img src="image_23.png" /><figcaption>**Figure 2.6:** Projection of the climate data volume increase, from Overpeck et al (2011).</figcaption></figure>

understand GAN is to see the Discriminator as the loss function of the Generator. Here again, various types of GAN exist, but all rely on the same principle.

## 2.6 Deep Learning and climate applications

The application of artificial neural networks to meteorological or climate problems is not new. Since the early 90s, several studies have used neural networks for forecasting weather or for statistical downscaling (Gardner and Dorling, 1998). Nevertheless, the recent boom of deep neural networks has also reached meteorology and climate science with an exponential increase in publications. Climate science has generated huge amounts of data from in-situ and satellite observations and climate models. The CMIP5 exercise generated about 5 Petabytes⁵ and CMIP6 about 23PB<fnref n="6" />, so about 5 times more in less than a decade. The projection made by Overpeck et al (2011) shown in Figure 2.6 seems to hold until now, and they forecast a stronger increase for the next decade. From a machine learning point of view, the strong advantage of climate data is that the climate modelling community produces clean, continuous and homogeneous gridded data. The outstanding performance of neural networks to deal with images notably but more generally to estimate non-linear relationships has brought much interest in the climate and weather communities. The convolutional layers are particularly adapted to climate maps since the different variables can be considered as the different channels of an image. In the following application examples, I will focus only on recent studies based on recent architecture.

Numerous studies have explored the capability of neural networks to forecast

<footnote n="51">⁵¹ PB=1 000 TB=1 000 000 GB</footnote>

<footnote n="6">⁶source: ESGF, data statistic page, Cinquini et al (2014)</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure 2.7:** Precipitation nowcasting Ravuri et al (2021b).</figcaption></figure>

weather at different time scales. Ravuri et al (2021a) recently presented a complex GAN architecture for precipitation nowcasting. The idea is to use live radar observation of precipitation and project the evolution of the precipitation objects for the following hours. The authors compare their method with a state-of-the-art stochastic approach (Pulkkinen et al, 2019) and two neural networks: one based on LSTM (Ayzel et al, 2020) and a second one based on a UNet model (Sønderby et al, 2020). The GAN network showed promising results and notably produced realistic maps with precise locations for precipitation, while the two other deep learning approaches create blurry images after some time steps. This result, illustrated in Figure 2.7, is typical of what GAN can offer: it is one possible outcome from the starting point while UNet, for example, tries to cover all the possibilities, which are increasing with time. Forecasting at other time scales has also been explored. Rasp et al (2020) proposed a benchmark dataset to evaluate and compare a data-driven medium-range (up to a week) forecasting model, including baseline scores from state-of-the-art numerical models and multiple studies already used it (Rasp and Thuerey, 2021; Scher and Messori, 2021). Among others, Lu et al (2022) or Martin et al (2022) have explored the predictability at sub-seasonal and seasonal scales, notably by focusing on mid-term oscillations. Recently, multiple review studies have tried to analyse the role of recent deep learning approaches in weather forecasting as it is not clear yet how far those methods will go (Schultz et al, 2021; Dewitte et al, 2021). If it seems evident that deep learning methods will help improve some aspects of the forecast, it is more difficult to say if an end-to-end approach will ever replace the current systems.

A second neural network application in the climate field is detecting and classifying meteorological objects. Indeed as mentioned earlier, climate fields can be considered images, and the good properties of convolutional neural networks in ob-
ject detection could help in many cases. For example, Racah et al (2017), introduce a labelled dataset from climate simulations and a baseline CNN architecture to encourage the development of an automatic and efficient method to detect in climate simulation extreme weather patterns. The labelled objects are Tropical Depression, Tropical Cyclones, Extra-Tropical Cyclones and Atmospheric Rivers, all subject to severe cause damage to human societies. Kim et al (2019) or Giffard-Roisin et al (2020) also propose NN-based method to track hurricanes and cyclones. Mounier et al (2022) uses very high resolution simulation to detect bow-echo<sup>7</sup> in forecast ensembles. Nevertheless, object detection necessitates a large volume of good quality labelled data which limits their development strongly.

In the introductory chapter, I briefly presented the role of parametrisations. They are an empirical representation of processes that occur at finer scales than the grid resolution. These parametrisations are based on expert knowledge and might lead to model biases (Medeiros and Stevens, 2011). Some studies have proposed using neural networks to learn the actual physical processes within high-resolution models that resolve those processes (Rasp et al, 2020; Brajard et al, 2021). The results are auspicious, and most research still has to be done as some stabilisation problems persist (Balogh et al, 2022).

Finally, as mentioned in the first chapter, neural network architecture has also been used in empirical statistical downscaling. Indeed, statistical downscaling aims to learn a generally non-linear relationship between GCM output (or similar large-scale representation) and station data or high-resolution gridded observational products. Convolutional neural networks have proven to deal with high dimensional data in input and output. Their ability to reconstruct images, as in super-resolution exercises, is also very interesting, as spatial downscaling aims to recreate high-resolution maps. Moreover, CNNs manage, by construction, to identify the right combination of inputs for each output (station or grid point). The use of the artificial neural network is not new in statistical downscaling (Cannon, 2008; Chadwick et al, 2011), and the recent performance of CNN has motivated many studies to use them in ESD. Baño-Medina et al (2020) have notably built simple CNN architecture to downscale precipitation and temperature over Europe in a perfect prognosis context with similar or better (depending on the metric) performance than state-of-the-art ESD regression. From the existing literature I see two ways of applying deep neural networks to ESD problems. First, in perfect prognosis approach as Baño-Medina et al (2020) or Höhlein et al (2020) following a classical regression task. Second, some studies have used them in a MOS approach, trying to link the low-resolution and the high-resolution from the same variables. The task is then similar to super-resolution problems that aim to increase the resolution of images. Wang et al (2021), Serifi et al (2021) and Babaousmail et al (2021) build RCM-emulators following this super resolution approach. In the following chapters, the proposed RCM-emulator is build in a perfect prognosis. Finally, it is remarkable that multiple studies have used a UNet architecture to downscale surface variables with generally promising results (Höhlein et al, 2020; Serifi et al, 2021; Sha et al, 2020a,b).

<footnote n="7">⁷Convective system responsible for extreme winds or tornados</footnote>
## 2.7 Link with this work

This second chapter was somehow an extension of the introductory chapter as it did not aim to present the results of my work but to gently introduce convolutional neural networks and justify why they are adapted for the RCM-emulator. In preliminary work for this PhD, I tried other machine learning methods. The aim was to explore the idea's feasibility, so I concentrated on single grid points. I used classical ESD methods such as linear regression or advanced quantile mapping (CDFt, Michelangeli et al, 2009) and various machine learning algorithms based on decision trees like Random Forest or XGboost. In all cases, the downscaling relationship seemed well estimated, which gave us confidence in the feasibility of the RCM-emulator. Nevertheless, the efficiency of those methods is questioned. Indeed they are all univariate in the sense that they can consider only a single output. Therefore, it is necessary to fit one model for each grid point which can get really heavy as we work with high-resolution models. Convolutional Neural Networks present the nice advantage of dealing efficiently with large images implying a single model for the emulator. Moreover, we saw in this chapter that CNN would ensure some spatial coherence as the convolution layers consider the image blockwise. Multiple studies have shown good results for CNN in similar exercises as the one I propose here. I compared CNN with the other methods over single grid points, and they showed similar or better results. Chapter 3 includes such a comparison over large domain. Furthermore, after some small tests, I decided to base my work on the UNet architecture. Chapter 3 gives more details about the specific architecture I built. However, it is important to keep in mind here that the network architecture does not guarantee anything. The construction of the training set, the choices regarding some parameters or the loss function also have a crucial influence. The following chapters will explore and discuss those choices.
# Chapter 3

## Regional Climate Model emulator based on deep learning: concept and first evaluation of a novel hybrid downscaling approach

A. Doury, S. Somot, S. Gadat, A. Ribes, L. Corre

This chapter is a published article in *Climate Dynamics*. I did not modify the corpus here. The organisation is the same, with the appendices and supplementary material at the end of the chapter. I simply modified the layout such that it stays coherent with the rest of the document. Some elements of the introduction and methodology section can be redundant with the previous chapters. The original article DOI is 10.1007/s00382-022-06343-9.
**3.1 Introduction**

Climate models are an essential tool to study possible evolutions of the climate
according to different scenarios of greenhouse gas emissions. These numerical mod-
els represent the physical and dynamical processes present in the atmosphere and
their interactions with other components of the Earth System. The complexity of
these models involves compromises between the computational costs, the horizontal
resolution and, in some cases, the domain size.

Global Climate Models (GCMs) produce simulations covering the whole planet at reasonable cost thanks to a low spatial resolution (from 50 to 300 km). The large number of different GCMs developed worldwide allows to build large and co-ordinated ensembles of simulations, thanks to a strong international cooperation. These big ensembles (CMIP3/5/6, Meehl et al, 2007; Taylor et al, 2012; Eyring et al, 2016) are necessary to correctly explore the different sources of variability and uncertainties in order to deliver reliable information about future climate change at large spatial scales. However, the resolution of these models is too coarse to derive any fine scale information, which is of primary importance for impact studies and adaptation policies. Consequently, it is crucial to downscale the GCM outputs to a higher resolution. Two families of downscaling have emerged: empirical-statistical downscaling and dynamical downscaling. Both approaches have their own strengths and weaknesses.

Empirical Statistical Downscaling methods (ESD) estimate functions to link large scale atmosphere fields with local scale variables using observational data. Local implications of future climate changes are then obtained by applying these functions to GCM outputs. Gutiérrez et al (2019) present an overview of ESD methods and evaluate their ability to downscale historical GCM simulations. The great advantage of these statistical methods is their computational efficiency, which makes the downscaling of large GCM ensembles possible. On the other hand, they have two main limitations due to their dependency on observational data. First of all, they are applicable only for regions and variables for which local long-term observations are available. Secondly, they rely on the stationary assumption of the large-scale / local-scale relationship, which implies that a statistical model calibrated in the past and present climate remains reliable in the forthcoming climate. Studies tend to show that the calibration period has a non-negligible impact on the results (Wilby et al, 1998; Schmith, 2008; Dayon et al, 2015; Erlandsen et al, 2020).

Dynamical Downscaling (DD) is based on Regional Climate Models (RCMs). These models have higher resolution than GCMs (1 to 50 km) but are restricted to a limited area domain to keep their computational costs affordable. They are nested in a GCM, e.g., they receive at small and regular time intervals dynamical information from this GCM at their lateral boundaries. One key advantage of RCMs is to rely on the same physical hypotheses as the one involved in GCMs. They provide a complete description of the state of the atmosphere over their domain through a large set of variables at high temporal and spatial resolution. The added value of RCMs has been demonstrated in several studies (Giorgi et al, 2016; Torma et al, 2015; Prein et al, 2016; Fantini et al, 2018; Kotlarski et al, 2015, for examples). In
order to deliver robust information about future local responses to climate change, it is necessary to explore the uncertainty associated with RCM simulations. Déqué et al (2007) and Evin et al (2021) assess that four sources of uncertainty are at play in a regional climate simulation: the choice of the driving GCM, the greenhouse gas scenario, the choice of the RCM itself and the internal variability. Their relative importance depends on the considered variables, spatial scale, and timeline. According to these results, it is important (Déqué et al, 2012; Evin et al, 2021; Fernández et al, 2019) to complete 4D matrices [SCENARIO, GCM, RCM, MEMBER] to deliver robust messages, where members are several simulations of each (SCENARIO, GCM, RCM) triplet. However, the main limitation of RCM is their high computational costs which makes the completion of such matrices impossible.

This study proposes a novel hybrid downscaling method to enlarge the size of RCM simulation ensembles. The idea is to combine the advantages of both dynamical and statistical downscaling to tackle their respective limits. This statistical RCM-emulator uses Machine Learning methods to learn the relationship between large scale fields and local scale variables inside regional climate simulations. It aims to estimate the downscaling function included in the RCM in order to apply it to new GCM simulations. This framework will allow to learn this function on the entire RCM domain, in past and future climate, under different scenarios. Besides, the emulator relies on Machine Learning algorithms with low computational costs, which will enable to increase RCM simulation ensembles and to better explore the uncertainties associated with these high resolution simulations.

Hybrid statistical-dynamical downscaling methods have already been proposed. They are methods which combine, in different ways, regional climate models and statistical approaches to obtain local climate information. Several studies such as Pryor and Barthelmie (2014), Vrac et al (2012) or Turco et al (2011) perform 2-step downscaling by applying ESD methods to RCM simulations. Colette et al (2012) apply bias correction methods to GCM outputs before downscaling with RCMs. Maraun and Widmann (2018) are among the first to mention the concept of emulators. Few studies have combined ESD and DD for the same purpose as in this study. For instance, Walton et al (2015) propose a statistical model which estimates from GCM outputs, the high resolution warming pattern for each month in California. It is calibrated using RCM simulations and relies on a simple linear combination of two predictors from the GCM (the monthly mean over the domain and an indicator for the land/sea contrast) plus an indicator for the spatial variance (obtained thanks to PCA). Berg et al (2015) adapt the same protocol for monthly changes in precipitation over California. With respect to those pioneer studies, we propose here to further develop this approach by using a neural network based method and by emulating the full time series at the input time scale, allowing to explore daily climate extremes.

In recent years, climate science has taken advantage of the recent strides in performances of Deep Learning algorithms (see Lecun et al, 2015, for an overview). Indeed, thanks to their capacity to deal with large amounts of data and the strong ability of Convolutional Neural Network (CNN) (LeCun et al, 1998) to extract high level features from images, these algorithms are particularly adapted to climate and
weather sciences. Reichstein et al (2019) present a complete overview of the climate studies applying Deep Learning and future possibilities. In particular, Vandal et al (2017, 2019) and Baño-Medina et al (2020, 2021) showed the good ability of Convolutional Neural Network (CNN) architecture to learn the transfer function between large scale fields and local variables in statistical downscaling applications. Baño-Medina et al (2021) confirms the suitability of CNN for ESD in an inter-comparison study, while Vandal et al (2017) demonstrates the good performances of CNN in front of a state-of-the-art bias correction method. The concept of emulator is mentioned in Reichstein et al (2019) as surrogate models trained to replace a complex and computationally expensive physical model (entirely or only parts of it). Once trained, this emulator should be able to produce simulations much faster than the original model. In this context, the RCM-emulator proposed here is based on a different fully convolutional neural network architecture known as UNet (Ronneberger et al, 2015). Wang et al (2021) propose a different emulator following a different strategy than ours as they train it using a low and high-resolution version of the same RCM and another type of neural network (namely CGAN).

This study presents and validates the concept of statistical RCM-emulator. We will focus on emulating the near-surface temperature in a RCM over a specific domain, including high mountains, shore areas, and islands in Western Europe. This domain regroups areas where the RCM presents added value compared to GCM but remains small enough to perform quick sensitivity tests. This paper is organized as follows: Section 2 presents the whole framework to define, train and evaluate the emulator, while Section 3 shows the emulator results. Finally, Sections 4 and 5 discuss the results of the emulator and provide conclusions.

## 3.2 Methodology

This section provides a complete description of the UNet-based RCM-emulator used in this paper. The notations are summarised in Table 3.1. The RCM-emulator uses a neural network architecture to learn the relationship between large-scale fields and local-scale variables inside regional climate simulations. RCMs include a downscaling function (F) which transforms large scale information (X,Z) into high resolution surface variables (Y). The statistical RCM-emulator aims to solve a common Machine Learning problem

$$Y = F(X, Z)$$

which is to estimate F by $\hat{F}$ in order to apply it to new GCM simulations. The following paragraphs describe the list of predictors used as inputs and their domain, the predictand (or target) and its domain, the neural network architecture, the framework used to train the emulator and the metrics used to evaluate its performances.

### 3.2.1 Models and simulations

This study focuses on the emulation of the daily near-surface temperature from EURO-CORDEX simulations based on the CNRM-ALADIN63 regional climate
<table><caption>**Table 3.2:** List of predictors</caption><thead><tr><th colspan="6">2D Variables</th></tr><tr><th>Field</th><th>Altitude Levels</th><th>Variables Notation</th><th>Units</th><th>Temporal Aggregation</th><th>Dimension</th></tr></thead><tbody><tr><td rowspan="3">Geopotential</td><td rowspan="3">850, 700, 500 hPa</td><td>zg500</td><td rowspan="3">m</td><td rowspan="3">Daily mean</td><td rowspan="3">[i, j]</td></tr><tr><td>zg700</td></tr><tr><td>zg850</td></tr><tr><td rowspan="3">Specific Humidity</td><td rowspan="3">850, 700, 500 hPa</td><td>hus500</td><td rowspan="3">K</td><td rowspan="3">Daily mean</td><td rowspan="3">[i, j]</td></tr><tr><td>hus700</td></tr><tr><td>hus850</td></tr><tr><td rowspan="3">Temperature</td><td rowspan="3">850, 700, 500 hPa</td><td>ta500</td><td rowspan="3">m/s</td><td rowspan="3">Daily mean</td><td rowspan="3">[i, j]</td></tr><tr><td>ta700</td></tr><tr><td>ta850</td></tr><tr><td rowspan="4">Eastward Wind</td><td rowspan="4">850, 700, 500 hPa</td><td>ua500</td><td rowspan="4">m/s</td><td rowspan="4">Daily mean</td><td rowspan="4">[i, j]</td></tr><tr><td>ua700</td></tr><tr><td>ua850</td></tr><tr><td>ua850</td></tr><tr><td rowspan="4">Northward Wind</td><td rowspan="4">850, 700, 500 hPa</td><td>va500</td><td rowspan="4">m/s</td><td rowspan="4">Daily mean</td><td rowspan="4">[i, j]</td></tr><tr><td>va700</td></tr><tr><td>va850</td></tr><tr><td>va850</td></tr><tr><td>Sea Level Pressure</td><td>Surface</td><td>psl</td><td>Pa</td><td>Daily mean</td><td>[i, j]</td></tr><tr><td>Total Aerosol Optical Depth forcing</td><td>Surface</td><td>TAOD</td><td></td><td>Monthly mean</td><td>[i, j]</td></tr></tbody><thead><tr><th colspan="6">1-D Variables</th></tr></thead><tbody><tr><td>Daily spatial means of 2D variables</td><td>x̄<br/>with x ∈ V2D</td><td></td><td></td><td>Daily</td><td>[#V2D]</td></tr><tr><td>Daily spatial standard deviation of 2D variables</td><td>x̄<br/>with x ∈ V2D</td><td></td><td></td><td>Daily</td><td>[#V2D]</td></tr><tr><td>Total Anthropogenic greenhouses gas forcings</td><td>ant_ghg</td><td></td><td></td><td>Yearly</td><td>[1]</td></tr><tr><td>Solar and Ozone forcings</td><td>sol, oz</td><td></td><td></td><td>Yearly</td><td>[2]</td></tr><tr><td>Seasonal Indicators<br/>Cos(2πt/365) ; Sin(2πt/365)</td><td>cos,sin</td><td></td><td></td><td>Daily</td><td>[2]</td></tr></tbody></table>

dictors. These variables are normalised (see Equation 3.1) according to their daily spatial mean and standard deviation so that they all have the same importance for the neural network before the training.

The set of 1D variables includes external forcing also given to the RCM: the total concentration of greenhouse gases and the solar and ozone forcings. It also includes a cosinus, sinus vector to encode the information about the day of the year. Given that the 2D variables are normalised at each time step by their spatial mean, they don't carry any temporal information. For this reason, the daily spatial means and standard deviations time series for each 2D variable are included in the 1D input, bringing the size of this vector to 43. In order to always give normalised inputs to the emulator, Z is normalized (see Equation 3.2) according to the means and standard deviations over a reference period (1971-2000 here) chosen during the emulator training. The same set of means and standard deviations will be used to normalise any low resolution data to be downscaled by the emulator.
$$ \tilde{X}_{t,i,j,x} = \frac{X_{t,i,j,x} - \bar{X}_{t,x}}{\dot{X}_{t,x}} \qquad (3.1) $$

$$ \text{with } \bar{X}_{t,x} = \frac{\sum_{i=1}^{I} \sum_{j=1}^{J} X_{t,i,j,x}}{I \times J} \quad \text{and} \quad \dot{X}_{t,x} = \sqrt{\frac{\sum_{i=1}^{I} \sum_{j=1}^{J} (X_{t,i,j,x} - \bar{X}_{t,x})^2}{I \times J}} $$

This decomposition of the large scale information consists in giving separately the spatial structure of the atmosphere (X) and the temporal information (Z) to the emulator. Thanks to the neural network architecture described in Section 3.2.4, we force the emulator to consider both items of information equally.

$$ \tilde{Z}_{t,z} = \frac{Z_{t,z} - \bar{Z}_{ref,z}}{\dot{Z}_{ref,z}} \qquad (3.2) $$

$$ \text{with } \bar{Z}_{ref,z} = \frac{\sum_{t \in R} Z_{t,z}}{nb \text{ days in } R} \quad \text{and} \quad \dot{Z}_{ref,z} = \sqrt{\frac{\sum_{t \in R} (Z_{t,z} - \bar{Z}_{ref,z})^2}{nb \text{ days in } R}} $$

where R is the reference period.

The aim of the emulator is to downscale GCM simulations. Klaver et al (2020) shows that the effective resolution of climate models is often larger (about 3 times) than its nominal resolution. For instance, CNRM-CM5 is more reliable at a coarser resolution, probably about 450 to 600 km, than at its own horizontal resolution (≈ 150 km). For this reason, the set of 2D predictors are smoothed with a 3 × 3 moving average filter. The grid of the GCM is conserved, with each point containing a smoother information than the raw model outputs.

For this study the input domain is defined around the target domain (described in Section 3.2.3). It is a 16 × 16 ($J = I = 16$ in Table 3.1) CNRM-CM5 grid box visible on Figure 3.1. Each observation given to the emulator (see Figure 3.1 for an illustration) is a day $t$ and it is composed of a 3D array ($X_{i,j,x}$), where the two first dimensions are the spatial coordinates and the third dimension lists the different variables chosen as predictors, and a 1D array ($Z_z$) regrouping all the 1 dimensional inputs.

### 3.2.3 Predictands

In this study, to assess the ability of the RCM-emulator to reproduce the RCM downscaling function, we focused on the emulation of the daily near surface temperature over a small but complex domain. The target domain for this study is a box of 64x64 RCM grid points at 12km resolution (about 600 000 km²) centred over the south of France (Figure 3.1). It gathers different areas of interest for the regional climate models. It includes three mountain massif (Pyrenees, Massif Central and French Alps) which are almost invisible at the GCM scale (specially the Pyrenees). The domain also includes coastlines on the Mediterranean side and on the Atlantic
<figure><img src="image_25.png" /><figcaption>**Figure 3.1:** Illustration of an observation for a randomly-chosen day. Left: each map represents a 2D input variables (X), on the input domain, and the blue numbers correspond to the 1D variables (Z). Right: an example of Y, the near surface temperature on the target domain.</figcaption></figure>

side. Thanks to a better representation of the coastline at the RCM resolution, it takes better into account the sea effect on the shore climate. It was also important for us to add small islands on our evaluation domain, such as the Baleares (Spain), since they are invisible on the GCM grid and the RCM brings important information. Finally, three major rivers (plotted in blue in Figure 3.1) are on the domain with interesting impacts on climate (commented in Section 3.3): the Ebro in Spain, the Garonne in southwest of France and the Rhone on the east of the domain. This domain should therefore illustrate the added value brought by a RCM at local scale and be a good test-bed on the feasibility of emulating high resolution models.

### 3.2.4 Deep learning with UNet

#### Neural network model as a black-box regression model

The problem of statistical downscaling and of emulation of daily near-surface variables may be seen as a statistical regression problem where we need to build the best relationship between the output response Y and the input variables (X, Z). When looking at the L² loss between the prediction Ŷ and the true Y, the optimal link (denoted by F below) is theoretically known as the conditional expectation:

$$F(X, Z) = E[Y | (X, Z)]. \qquad (3.3)$$
Unfortunately, since we will only have access to a limited amount of observations collected over a finite number of days, we shall work with a training set formed by the collected observations ($(X_t, Z_t), Y_t)_{1\le t\le T}$ and try to build from the data an empirical estimation $\hat{F}$ of the unknown $F$.

For this purpose, we consider a family of relationship between $(X, Z)$ and $Y$ generated by a parametric deep neural network, whose architecture and main parameters are described later on. We use the symbol $\theta$ to refer to the values that describe the mechanism of one deep neural network, and $\Theta$ as the set of all possible values of $\theta$. Hence, the family of possible relationships described by a collection of neural networks correspond to a large set $(F_\theta)_{\theta\in\Theta}$. Deep learning then corresponds to the minimization of a $L^2$ data-fidelity term associated to the collected observations:

$$ \hat{F} = \arg \min_{F_\theta, \theta \in \Theta} \sum_{t=1}^{T} \|Y_t - F_\theta(X_t, Z_t)\|^2. \qquad (3.4) $$

## Training a deep neural network

To train our emulator between low resolution fields and one high resolution target variable, we used a neural network architecture called UNet whose architecture is described below. As usual in neural networks, the neurons of UNet are organised in layers. Given a set $E_n$ of input variables denoted by $(x_i)_{i\in E_n}$ of an individual neuron $n$, the output of $n$ corresponds to a non-linear transformation $\phi$ (called activation function) of a weighted sum of its inputs:

$$ \mathbf{n}((x_i)_{i \in E_n}) = \phi(\langle\langle(w_i)_{i \in E_n}, (x_i)_{i \in E_n}\rangle\rangle) = \phi\left(\sum_{i \in E_n} w_i x_i\right). \qquad (3.5) $$

The connection between the different layers and their neurons then depends on the architecture of the network. In a fully connected network (multilayer perceptron) all the neurons of a hidden layer are connected to all the neurons of the previous layer. The deepness of a network then depends on the number of layers.

As indicated in the previous paragraph, the machine learning procedure corresponds to the choice of a particular set of weights over each neuron to optimise a data fidelity term. Given a training set, a deep learning algorithm then solves a difficult multimodal minimisation problem as the one stated in (3.4) with the help of gradient descent strategies with stochastic algorithms. The weights associated with each neuron and each connection are then re-evaluated according to the evolution of the loss function, following the backpropagation algorithm Rumelhart and Hinton (1986). This operation is repeated over all the examples until the error is stabilized. Once the neural network is trained, it may be used for prediction, i.e., to infer the value $Y$ from new inputs $(X, Z)$.

We emphasize that the bigger the dimension of the inputs and outputs, the larger the number of the parameters to be estimated and so the bigger the training set must be. Therefore, the quality of the training set is crucial: missing or wrong values will generate some additional fluctuation and errors in the training process. Moreover, we also need to cover a sufficiently large variety of scenarios in the input variables to ensure that our training set covers a wide range of possible inputs. For all these reasons, climate simulation datasets are ideal to train deep learning networks.
<figure><img src="image_26.png" /><figcaption>**Figure 3.2:** Scheme of the neural network architecture used for UNet emulator. The part of the network in the red frame corresponds to the original UNet defined in Ronneberger et al (2015)</figcaption></figure>

### UNet architecture

The emulator proposed in this study relies on a specific architecture called UNet introduced by Ronneberger et al (2015). It is a fully convolutional architecture, i.e. all layers are convolutional layers known for their strength to deal efficiently with image data. UNet is known for its good ability to identify different objects and areas in an image. It involves gathering pixels that correspond to the same object. This key feature is naturally interesting for meteorological maps. The emulator needs to identify the different meteorological structures present in the low resolution predictors for a given day in order to predict the corresponding high resolution near-surface temperature. Moreover, the good performances of UNet have been demonstrated in problems such as pixel-wise regression (Yao et al, 2018) or super-resolution (Hu et al, 2019), which are comparable to ours.

The original UNet architecture is illustrated in the red frame in Figure 3.2 and more precisely described in appendix 3.6.2. It is composed of an encoding part on the left side, which reduces the spatial dimension of the input image while the pixel dimension gets deeper. On the right side, the expansive path decodes the encoded information to reconstruct the target image. Each encoding step is recalled in the decoding part, allowing the network to find the best way to the target image through different levels.

The architecture we propose is adapted from the original UNet to suit our problem. Firstly, the input object is a description of the climate conditions leading to the target image. It is then a 3D object while the classical UNet takes a single image as input. This aspect differs from the natural use of UNet but does not modify its construction. Secondly, we added a second source of input, the set of 1D predictors (see table 3.2) which is given at the bottom of the UNet (see Figure 3.2). This 1D
<figure><img src="image_27.png" /><figcaption>**Figure 3.3:** Time series of spatial correlation of the atmospheric temperature at 700hpa between ALADIN63 and its driving GCM, CNRM-CM5, over the input domain.</figcaption></figure>

object is encoded by a dense network and then concatenated with the encoded 2D inputs. The "U" shape of the UNet architecture allowed us to give this information such that the network treats equally the two sources of input. Finally, we extended the expansive part of the UNet to reach the higher resolution of the output as illustrated in Figure 3.2.

We chose the Mean Square Error (MSE) as the loss function to train the network, as we have a regression problem. Moreover, it is well adapted for variables following Gaussian-like distribution, such as temperature. The neural network was built and trained using the Keras Tensorflow environment in Python (https://keras.io). The network trained for this study has about 25 million parameters to fit.

### 3.2.5 Training of the emulator : Perfect Model Framework

As any statistical downscaling and any machine learning method, the emulator needs to be calibrated on a training set. It consists in showing the emulator many examples of predictors and their corresponding target such that the parameters of the network can be fitted as mentioned in Section 3.2.4. The emulator is trained in a perfect model framework, with both predictors and predictands coming from the same RCM simulation. The intuitive path to train the emulator is to use GCM outputs as predictors and its driven RCM outputs as predictands, but there are many reasons for our choice. First of all, it guarantees a perfect temporal correlation between large scale predictors and a local scale predictand. Indeed, Figure 3.3 shows that GCM and RCM large scales are not always well correlated, with an average correlation of 0.9 and 10% of the days with a coefficient of correlation lower than 0.75. These mismatches are quite well known and often due to internal variability as explained by Sanchez-Gomez et al (2009); Sanchez-Gomez and Somot (2018). Moreover, there are more consistent biases (discussed in Section 3.4.1) between GCM and RCM large scales. It is of primary importance that the inputs and outputs used to calibrate the model are perfectly consistent, otherwise, the emulator will try to learn a non-existing or non-exact relationship. In this context, the perfect model framework allows us to focus on the downscaling function of the RCM, specifically. This approach is similar to the "super-resolution downscaling" mentioned by Wang et al (2021) and deployed by Vandal et al (2017) using observational data in an empirical statistical downscaling framework.

The training protocol is summarised in Figure 3.4. In a first step, the RCM
<figure><img src="image_28.png" /><figcaption>Figure 3.4: Scheme of the protocols for the training (left) and the two steps of evaluations (center and right).</figcaption></figure>

simulation outputs are upscaled to the GCM resolution (about 150km) thanks to a conservative interpolation. This first step transforms the RCM outputs into GCM-like outputs. This upscaled RCM is called UPRCM in the rest of this paper. In the second step, these UPRCM outputs are smoothed by a 3x3 moving average filter to respect the protocol described in Section 3.2.2. This smoothing also targets to delete any local scale information which might persist through the upscaling step (as discussed in Section 3.4.3).

The near-surface temperature on the target domain is extracted from the same RCM simulation. Following this procedure, the emulator is trained using the ALADIN63 simulation forced by CNRM-CM5, covering the 1950-2100 period with the RCP8.5 scenario from 2006. As shown in Figure 3.4, we will use a different simulation driven by a different RCP scenario for the evaluation. We chose the two most extreme simulations (historical and RCP8.5) for the training in order to most effectively cover the range of possible climate states since the emulator does not target to extrapolate any information. Future studies could explore the best combination of simulations to calibrate the emulator.

### 3.2.6 Evaluation metrics

The emulated temperature series ($\hat{Y}$) will be compared to the original RCM series (Y) (mentioned as “RCM truth” in the rest of this paper) through statistical scores described below. Each of these metrics will be computed in each point over the complete series:

* **RMSE.** The Root Mean Squared Error measures the standard deviation of the prediction error (in °C):

$$RMSE(Y, \hat{Y}) = \sqrt{\frac{1}{T} \sum_{t} (Y_t - \hat{Y}_t)^2} \qquad (3.6)$$

* **Temporal Anomalies Correlation.** This is the Pearson correlation coefficient after removing the seasonal cycle:

$$ACC(Y, \hat{Y}) = \rho(Y_a, \hat{Y}_a), \qquad (3.7)$$
with $\rho$ the Pearson correlation coefficient and $Y_a$ and $\hat{Y}_a$ are the anomaly series after removing a seasonal cycle computed on the whole series.

* **Ratio of Variance.** It indicates the performance of the emulator in reproducing the local daily variability. We provide this score as a percentage.

$$ RoV(Y, \hat{Y}) = \frac{Var(\hat{Y})}{Var(Y)} * 100 \qquad (3.8) $$

* **Wasserstein distance.** It measures the distance between two probability density functions ($P, Q$). It relies on the optimal transport theory (Villani, 2009) and measures the minimum required “energy” to transform $P$ into $Q$. The energy here corresponds to the amount of distribution weight that is needed to be moved multiplied by the distance it has to be moved. In this study we use the 1-d Wasserstein distance, and its formulation between two samples becomes a rather simple function of ordered statistics:

$$ W_1(f(Y), f(\hat{Y})) = \sum_{i=1}^{T} |Y_{(i)} - \hat{Y}_{(i)}|, \qquad (3.9) $$

with $f(\bullet)$ the probability density function associated with the sample $\bullet$.

* **Climatology.** We compare the climatology maps over present (2006-2025) and future (2081-2100, not shown in Section 3.3) climate. The RCM truth and emulator maps are shown with their spatial correlation and RMSE. The error (emulator minus RCM) map is also computed.

* **Number of days over 30°C.** Same as climatology for the maps showing the number of days over 30°C.

* **99th Percentile.** Same as climatology for the maps showing the 99th percentile of the daily distribution.

* **Climate Change.** Climate change maps for the climatology, the number of days over 30°C and the 99th percentile (delta between future (2080-2100) and present (2006-2100) period).

These metrics are at the grid point scale and are presented as maps. However, to summarise these maps with few numbers we can compute their means and their super-quantile of order 0.05 (SQ05) and 0.95 (SQ95). The super-quantile $\alpha$ is defined as the mean of all the values larger (resp. smaller) than the quantile of order $\alpha$, when $\alpha$ is larger (resp. smaller) than 0.5. These values are shown in the Figures of Section 3.3 and Tables 3.3 and 3.4 in supplementary material.

### 3.2.7 Benchmark

For this study, we propose as benchmark the near surface temperature from the input simulation (before the moving average smoothing), interpolated on the target grid by bilinear interpolation. As this study is the first to propose such an emulator, there is no already established benchmark. The one proposed here is a naive high-resolution prediction given available predictors (low-resolution), it is used as a basic
reference and not a potential competitor for the emulator. It allows the reader to
locate the emulator somewhere in between the simplest possible downscaling (simple
interpolation of the original low resolution simulation) and the most complex one
(RCM simulation). All the metrics introduced in Section 3.2.6 will be applied to
our benchmark.

## 3.3 Results

This section presents the emulator performances in terms of its computational
costs (in Section 3.3.1) and its ability to reproduce the near surface temperature time
series at high resolution. As illustrated in Figure 3.4, we will evaluate the emulator in
two steps, (1) in the perfect model world (Section 3.3.2) and (2) when the emulator
inputs come from a GCM simulation (Section 3.3.3). The RCM simulation used to
evaluate the model (also called target simulation) is the ALADIN63, RCP45, 2006-
2100 forced by the GCM simulation CNRM-CM5, RCP45, 2006-2100. Note that the
emulator never saw the target simulation during the training phase. This evaluation
exercise illustrates a potential use of the emulator: downscaling a scenario that the
RCM has not previously downscaled.

### 3.3.1 Computational efficiency

The emulator is trained on a GPU (Nvidia, GeForce GTX 1080 Ti). About 60
epochs are necessary to train the network, and each epoch takes about 130 seconds
with a batch size of 100 observations. The training of the emulator takes about
2 hours. Once the emulator is trained, the downscaling of a new low resolution
simulation is almost instantaneous (less than a minute). It is a significant gain in
time compared to RCM, even if these time lengths do not include the preparation of
the inputs, which depends mainly on the size of the series to downscale and on the
data accessibility. It would, when including the input preparation, take only a few
hours on a simple CPU or GPU to produce a simulation with the trained emulator,
while it takes several weeks to perform a RCM simulation on a super-computer.

### 3.3.2 Evaluation step 1 : Perfect model world

In a first step, the emulator is evaluated in the perfect model world, meaning
that the inputs come from the UPRCM simulation. This first evaluation step is
necessary to control the performances of the emulator in similar conditions as during
its training. Moreover, the perfect model framework guarantees perfect temporal
matches between the large scale low resolution fields (the inputs) and the local scale
high resolution temperature (the target). The emulator should then be able to
reproduce perfectly the temperature series that is simulated by the RCM. This first
evaluation of the emulator is divided in two parts. In the section 3.3.2 we analyse the
ability of the emulator to reproduce the RCM simulation. The section 3.7 compares
the specific emulator proposed in this study with respect to two others emulators
relying on standard empirical statistical downscaling methods.
<figure><img src="image_29.png" /><figcaption>**Figure 3.5:** Randomly chosen illustration of the production of the emulator (in evaluation step 1) with inputs coming from UPRCM: (a) temperature (°C) at a random day over the target domain for the raw UPRCM, the interpolated UPRCM, the emulator and the RCM truth and (b) random year time series (°C) for 4 particular grid points.</figcaption></figure>

## Ability to reproduce the RCM simulation

The emulator aims to learn and reproduce the downscaling function included in the RCM, i.e., to transform the low resolution daily information about the state of the atmosphere into high-resolution daily surface temperature. We dedicate this first evaluation section to comparing the prediction of the emulator with the RCM truth in perfect model. The benchmark for this first evaluation is the UPRCM near surface temperature re-interpolated on the RCM grid. It is referred as “I-UPRCM”. We are aware that the “I-UPRCM” field is a very simple benchmark and can not compete with the emulator in terms of realism but it allows to simply measure the action of the emulator regarding the input simulation.

Figure 3.5a illustrates the production of the emulator for a random day regarding the target and the benchmark. The RCM truth map presents a refined and complex spatial structure largely missing in the UPRCM map. Moreover, it is evident on
<figure><img src="image_30.png" /><figcaption>**Figure 3.6:** (a) Daily probability density functions from the RCM truth, the emulator (in evaluation step 1) and the I-UPRCM at 4 particular grid points over the whole simulation period. (b) (resp. c)) Maps of performance scores of the emulator (resp. of the I-UPRCM) with respect to the RCM truth computed over the whole simulation period. For each map, the values of the spatial mean and super-quantiles (SQ05 and SQ95) are added.</figcaption></figure>
the I-UPRCM map that the simple bilinear interpolation does not recreate these
high resolution spatial patterns. The emulator shows for this given day an excellent
ability to reproduce the spatial structure of the RCM truth. It has very accurate
spatial correlation and RMSE and estimates the right temperature range. On Fig-
ure 3.5b, we show the daily time series for four specific points shown on the RCM
truth map (Marseille, Toulouse, a high Pyrenees grid point and a point in Majorca)
during a random year. The RCM transforms the large scale temperature (visible
on I-UPRCM) differently over the four points. In the Pyrenees, the RCM shifts
the series and seems to increase the variance. In Marseille, it appears to produce
a warmer summer without strongly impacting winter characteristics, and it seems
to be the exact opposite in Majorca. On the contrary, in Toulouse, I-UPRCM and
RCM are close. For each of these 4 cases, the emulator reproduces the sophistication
of the RCM series almost perfectly.

Figure 3.5a gives the impression that the emulator has a good ability to reproduce
the complex spatial structure brought by the RCM, and we can generalise this
result with the other Figures. First of all, the performance scores (Figure 3.6b)
of the emulator confirm that this good representation of the temperature’s spatial
structure is robust over the whole series. The spatial homogeneity of the score maps
tends to show that the emulator does not have particular difficulties over complex
areas. The spatial correlation (equal to 1) and the very low spatial RMSE (0.07°C)
of the climatology maps in Figure 3.7 in the present climate support this result. In
particular, it is worth noting on the climatology maps that the altitude dependency
is well reproduced as well as the warmer patterns in the Ebro and Rhone valley
or along the coastlines. Moreover, the comparison with the interpolated UPRCM
shows the added-value of the emulator and in particular its ability to reinvent the
fine-scale spatial pattern of the RCM truth from the large-scale field. Indeed, the
score maps of the I-UPRCM (Figure 3.6c) or the present climatology error map
( Figure 3.7a) shows strong spatial structures, highlighting the regions where the
RCM brings added value and that the emulator reproduces successfully.

The RCM resolution also allows to have a better representation of the daily vari-
ability at the local scales over critical regions. The difference in variance between the
I-UPRCM and the RCM is visible on Figure 3.6c. The I-UPRCM underestimates
the variability and is poorly correlated over the higher reliefs, the coastlines and
the river valley. In contrast, the emulator reproduces more than 90% of the RCM
variance over the whole domain. The RMSE and temporal correlation maps of the
emulator confirm the impression given by Figure 3.5b that it sticks almost perfectly
to the RCM truth series. Moreover, the RCM daily variability is strongly dependent
on the region. Indeed, the RCM transforms the “I-UPRCM” pdfs in different ways
across the domain (visible on Figures 3.6ac). Figures 3.6ab show that the emulator
succeeds particularly well in filling these gaps.
<figure><img src="image_31.png" /></figure>
**Figure 3.7:** (a) Maps of long-term mean climatologies, (b) number of days over 30°C and (c) 99th percentile of daily near-surface temperature for a present-climate period (2006-2025) and for the climate change signal (2080-2100 minus 2006-2025), for the RCM truth, the emulator (in evaluation step 1) and the I-UPRCM. On each line, the two last maps show the error map of the emulator and the I-UPRCM. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added, as well as the spatial correlation and spatial RMSE for the emulator and I-UPRCM maps.

The emulator’s good representation of the daily variability and temporal correlation involves a good representation of the extreme values. The probability density functions of the four specific points on Figure 3.6a show that the entire pdfs are fully recreated, including the tails. The Wasserstein distance map extends this result to the whole target domain. The two extreme scores computed for the present climate on Figure 3.7bc confirm these results. The 99th percentile emulated map is almost identical to the target one verified by the difference map, with a maximum difference of less than 1°C for values over 35°C. The spatial pattern of the 99th percentile map is here again correctly captured by the emulator, particularly along the Garonne river that concentrates high extremes. The number of days over 30°C is a relatively more complicated score to reproduce since it involves an arbitrary threshold. The emulator keeps performing well with a high spatial correlation between the emulator and the RCM truth. However, it appears that the emulator misses some extreme days, involving a lack in the intensity of some extremes metrics.

Finally, the high-resolution RCM produces relevant small-scale structures in the climate change maps. In particular, RCMs simulate an elevation-dependent warming (see the Pyrenées and Alps areas Kotlarski et al, 2015), a weaker warming near the coasts (see the Spanish or Atlantic coast) and a specific signal over the islands as shown in the second lines of Figures 3.7ac. It can be asked if the emulator can reproduce these local specificities for the climate change signal. The emulator is able to capture this spatial structure of the warming but with a slight lack of intensity which is general over the whole domain. The reproduction of the climate change in the extremes suffers the same underestimation of the warming but also offers the same good ability to reproduce the spatial structure, with high spatial correlation.

This first evaluation step shows that if the emulator is still perfectible, in particular when looking at extremes or climate change intensity, it is able to almost perfectly reproduce the spatial structure and daily variability of the near surface temperature in the perfect model world.

## Comparison to simpler emulators

This section validates the UNet-based emulator by comparing its performances to two other emulators based on standard empirical statistical downscaling methods. The perfect model framework allows comparing different emulators since the RCM truth is the ideal reference when downscaling the UPRCM simulation.
<figure><img src="image_32.png" /><figcaption>**Figure 3.8:** (a) (resp. b) Maps of performance scores of the CDFt emulator (resp. of the MLR emulator) with respect to the RCM truth computed over the whole simulation period in perfect model evaluation. For each map, the values of the spatial mean and super-quantiles (SQ05 and SQ95) are added. (c) Climate change maps and the difference maps with respect to the RCM truth for the CDFt emulator and the MLR emulator. As on figure 3.6 and 3.7 the spatial scores are added on the maps.</figcaption></figure>

The two comparison emulators proposed here rely on very classical ESD methods. The first one is CDFt (Michelangeli et al, 2009; Vrac et al, 2012), which belongs to the Model Output Statistic family of methods. It can be considered as an extended "quantile-quantile" method. It transforms the low resolution temperature cumulative distribution function into the high resolution one. The second method is a simple Multiple Linear Regression (MLR Huth, 2002; Huth et al, 2015), which belongs to the Perfect Prognosis family of method. The inputs used for the MLR are the same as for the UNet emulator, but only the closest low resolution
point is used for the 2D variables. The methods and the way we used them as RCM-emulator are more precisely described in appendix 3.6.1. Both methods are evaluated in the VALUE project (Gutiérrez et al, 2019), showing reasonable results. They are trained in the same conditions as the UNet emulator, following the perfect model framework described in Section 3.2.5.

The results are illustrated on Figure 3.8 and a more complete evaluation can be found in the appendix 3.6.1. As for the neural network emulator, both CDFt and MLR emulators present a good ability to reproduce the small scale information brought by the RCM. They present a good variance ratio and temporal correlation leading to a good RMSE. The climate change maps (Figure 3.8c) show that both methods are also able to reproduce the RCM climate change maps, even if the MLR method seems to have more difficulties and particularly over the relief. These results are expected as these methods showed good performances in statistical downscaling and both methods could constitute reasonable emulators. However, the neural network emulator performs better in most metrics. Indeed, the UNet Emulator presents a lower RMSE showing a better accuracy to fit the original RCM series. It also better captures the spatial structure of the RCM truth, as shown by the spatial correlation scores on the climatology, and extremes maps (Figures 3.7, 3.8 and appendix 3.6.1).

According to this comparison, we conclude that the three approaches may con-
stitute good RCM-emulators but the one based on the UNet approach is far
better, especially for the correlation with the RCM truth. As the UNet-approach
is not significantly more complex to apply and not costly to train, we decide
to select it as our main emulator for the following of the article, that is to say
to illustrate the way RCM-emulators can be applied in practical to downscale
GCM simulations. Furthermore, the UNet based emulator, as well as the two
simpler emulators, are possible emulator among many other options. We invite the
statistical downscaling community to propose other emulators, relying on different
statistical methods, trained in different ways, to find the best tool to emulate RCMs.

3.3.3 Evaluation step 2 : GCM world

In this second evaluation step, we directly downscale a GCM simulation. The benchmark for this evaluation is the near surface temperature from the GCM, interpolated on the target grid. It will be referred to as I-GCM. The emulated series and the benchmark are compared to the RCM simulation driven by the same GCM simulation.

Figure 3.9a illustrates the production of the emulator regarding the benchmark and RCM truth for the same day as Figure 3.5a. First of all, as for the I-UPRCM, the I-GCM map does not show any of the complex RCM spatial structures. The I-GCM is less correlated with the RCM and warmer than the I-UPRCM. In contrast, the emulator reproduces the complex spatial structure of the RCM very well with a spatial correlation of 0.98 but appears to have a warm bias with respect to the RCM truth. The four time series are consistent with the previous section, with
<figure><img src="image_33.png" /><figcaption>**Figure 3.9:** Randomly chosen illustration of the production of the emulator (in evaluation step 2) with inputs coming from the GCM: (a) near-surface temperature (°C) at a random day over the target domain for the GCM, the interpolated GCM, the emulator, and RCM truth and (b) random year time series (°C) for 4 particular grid points.</figcaption></figure>

fundamental differences between I-GCM and RCM, which the emulator captured very well. However, the correlation between the emulator and the RCM seems to be not as good as in the perfect model framework.

Figure 3.9b is a very good illustration of the RCM-GCM large scale de-correlation issue presented in Section 3.2.4. Indeed the less good correlation of the emulator with the RCM is probably due to mismatches between GCM and RCM large scales. For instance, in the beginning of November, on the time series shown on Figure 3.9b, the RCM seems to simulate a cold extreme on the whole domain, which appears neither in the interpolated GCM nor in the emulator. The same kind of phenomenon occurs regularly along the series and is confirmed by lower temporal correlations between the RCM truth and the I-GCM (Figure 3.10c) than with the I-UPRCM (Figure 3.6c). According to this, the emulated series can not present a good temporal correlation with the RCM truth since it is a daily downscaling of the GCM large scale. Keeping in mind these inconsistencies, it is still possible to
analyse the performances of the emulator if we leave aside these scores which are influenced by the poor temporal correlation (RMSE, ACC).

As in the first step of the evaluation (Section 3.3.2), the spatial structure of the RCM truth is well reproduced by the emulator. The present climatology map (Figure 3.11a) has a perfect spatial correlation with the RCM. The added value from the emulator is clear if compared to the interpolated GCM. The spatial temperature gradients simulated by the GCM seem to be mainly driven by the distance to the sea. On the other hand, the emulator manages to recreate the

<figure><img src="image_34.png" /><figcaption>**Figure 3.10:** (a) Daily probability density functions from the RCM truth, the emulator (in evaluation step 2) and the I-GCM at 4 particular grid points over the whole simulation period. (b) (resp. (c)) Maps of performance scores of the emulator (resp. of the I-GCM) with respect to the RCM truth computed over the whole simulation period. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added.</figcaption></figure>
complex structures created by the RCM, related to relief and coastlines. The emulator capacity to reproduce the RCM spatial structure seems as good as in Section 3.3.2. The scores in Figure 3.10b present a good spatial homogeneity, exactly like in the previous section 3.6b.

The error map in Figure 3.11(a) shows that the emulator is warmer in the present climate than the RCM truth (+0.96°C). This bias presents a North-South gradient with greater differences over the North of the target domain, which is consistent with the Wasserstein distance map on Figure 3.10b. The Wasserstein metric shows that the density probability functions from the emulated series are further away from the RCM truth with GCM-inputs than in UPRCM mode. The similarities between the Wasserstein scores and the present climatology difference map indicate that the emulator shifts the mean.

The daily variability is well reproduced by the emulator. As mentioned before, the weaker RMSE (Figure 3.10b) is mainly due to the lower correlation between GCM and RCM. But the ratio of variance demonstrates that the emulator manages to reproduce the daily variability over the whole domain. The RCM brings a complex structure of this variability (higher variability in the mountains than in plains, for example), and the emulator, as in the first evaluation step, recreates this fine scale. Moreover, the daily pdfs of the emulator (Figure 3.10a) are very consistent with the RCM ones, and the same range of values is covered for each of the four particular points.

This good representation of daily variability tends to suggest that the emu-
lator can reproduce the local extremes. Figures 3.11bc confirm these results,
with a very high spatial correlation between the emulator and the RCM truth
maps in present climate. The warmer extremes along the three river valleys
are present in both RCM and emulator maps, while they are absent from the
I-GCM maps. The warm bias observed in the present climatology map also
impacts these scores. The emulator map of the number of days over 30°C in
the present climate shows more hot days than the RCM but the same spatial
structure. The map of the 99th percentile over the 2006-2025 period shows the
same observation, with a warm bias (+0.82) slightly lower than the climatology bias.

Finally, the climate change signal is also well captured by the emulator. The different spatial patterns that bring the high resolution of the RCM in the Figures 3.11abc are also visible in the emulator climate change Figures. The emulator represents a weaker warming than the RCM, observable in average warming but also on the map of extremes. This underestimated warming is mainly due to the warm bias between GCM and RCM, which is less intense in the future. For instance, the warming from the emulator is 0.27°C weaker on average over the domain (with almost no spatial variation) than in the RCM. This number corresponds approximately to the cold bias from the I-GCM (0.19°C) plus the missed warming by the emulator in the perfect model framework (0.07). This tends to show that the emulator performs well in the GCM world but reproduces the GCM-RCM biases.
<figure><img src="image_35.png" /></figure>
**Figure 3.11:** (a) Maps of long-term mean climatologies, (b) number of days over 30°C and (c) 99th percentile of daily near-surface temperature for a present-climate period (2006-2025) and for the climate change signal (2080-2100 minus 2006-2025), for the RCM truth, the emulator (in evaluation step 2) and the I-GCM. On each line, the two last maps show the error map of the emulator and the I-GCM. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added, as well as the spatial correlation and spatial RMSE for the emulator and I-GCM maps.

This section shows that the emulator remains robust when applied to GCM inputs since it provides a realistic high-resolution simulation. As in the first step, the emulator exhibits several desirable features with an outstanding ability to reproduce the complex spatial structure of the daily variability and climatology of the RCM. We also showed that the emulator remains consistent with its driving large scale, which leads to inconsistencies with the RCM. In the next section, we will develop this discussion further.

# 3.4 Discussion

## 3.4.1 On the inconsistencies between GCM and RCM

Several recent studies (Sørland et al, 2018; Bartók et al, 2017; Boé et al, 2020) have highlighted the existence of large scale biases for various variables between RCMs and their driving GCM, and have discussed the reasons behind these inconsistencies. From a theoretical point of view, it is still controversial as to whether these inconsistencies are for good or bad reasons (Laprise et al, 2008) and therefore if the emulator should or should not reproduce them.

In our study, the emulator is trained in such a way that it focuses only on learning the downscaling function of the RCM, i.e., from the RCM large scale to the RCM small scale. Within this learning framework, the emulator can not learn GCM-RCM large-scale inconsistencies, if there should be any. Therefore, when GCM inputs are given to the emulator, the estimated RCM downscaling function is applied to the GCM large scales fields, and any GCM-RCM bias is conserved between the emulated series and the RCM one. Figure 3.12 shows the biases for the present-climate climatology between the GCM and the UPRCM over the input domain for TA700 and ZG700, at the GCM resolution. The GCM seems generally warmer than the UPRCM, which could partly explain the warm bias observed between the emulator results and the RCM truth in present climate (e.g., Figure 3.11).

These large scale biases between GCM and RCM raise the question of using the RCM to evaluate the emulator when applied to GCM data. Indeed, if these inconsistencies are for bad reasons (e.g., inconsistent atmospheric physics or inconsistent forcings), the emulator somehow corrects the GCM-RCM bias for the emulated variable. In this case, the RCM simulation cannot be considered as the targeted truth. However, if the RCM revises the large scale signal for good reasons (e.g., upscaling of the local added-value due to refined representation of physical processes), then the design of the emulator should probably be adapted.
<figure><img src="image_36.png" /><figcaption>**Figure 3.12:** Present (2006-2025) climatology differences for the atmospheric temperature and geopotential at 700 hpa: CNRM-CM5 RCP45 minus ALADIN63 driven by CNRM-CM5 RCP45 upscaled on the GCM grid.</figcaption></figure>

In future studies, we plan to use RCM runs with spectral nudging (Colin et al, 2010), two-way nested GCM-RCM runs or global high-resolution simulations for testing other modelling frameworks to further develop and evaluate the emulator.

### 3.4.2 On the stationary assumption

In the introduction of this paper, we state that the stationary assumption is one of the main limitations of empirical statistical downscaling. The emulator proposed here is similar in many ways to a classical ESD method, the main difference being that the downscaling function is learnt in a RCM simulation. The framework used to train the emulator is a good opportunity to test the stationary assumption for the RCM-emulator. We train the same emulator, with the same neural network architecture and same predictor set, but on the historical period (1951-2005) only. Results are reported in Tables 3.3 and 3.4 in the supplementary material, this version being named 'Emul-Hist'.

The perfect model (Table 3.3) evaluation constitutes the best way to evaluate the validity of this assumption properly. Emul-Hist has a cold bias over the whole simulation regarding the RCM truth and the range of emulators described in subsection 3.4.4. Moreover, this bias is much stronger for the future period (from 0.3°C in 2006-2025 to 1.3°C in 2080-2100). Emul-Hist manages to reproduce only 30% of the climate change simulated by the RCM. It also fails to capture most of the spatial structure of the warming since the spatial correlation between the Emul-hist and RCM climate change maps (0.86) is closer to the I-UPRCM (0.82) and than the main emulator (0.95) (see Figure 3.13). The Emul-hist average RMSE (1.35°C) over the whole series is also out of emulator range ([0.8; 0.86]). Results in GCM evaluation are also presented (Table 3.4), but due to the lack of proper reference, it
<figure><img src="image_37.png" /><figcaption>**Figure 3.13:** Maps of climatologies over present period (2006-2025) and climate change signal (2080-2100 versus 2006-2025), for the RCM truth, the Emulator presented in section 3.3 and the Emul-Hist. The two last columns correspond to the error maps of the Emulator and Emul-Hist with respect to the RCM truth. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added, as well as the spatial correlation and spatial RMSE for the Emul-Hist and I-UPRCM maps.</figcaption></figure>

is difficult to use them to assess the stationary assumption. However, it presents the same cold bias regarding the ensemble of emulators. These results demonstrate the importance of training the emulator in the wider range of possible climate states. We underline that not all ESD methods are expected to behave that poorly with respect to projected warming. However learning in the future is one of the main differences between our RCM emulator approach and the standard ESD approach that relies on past observations.

### 3.4.3 On the selection of the predictors

For this study, we chose to use a large number of inputs with almost no prior selection, leaving the emulator to select the right combination of inputs for each grid point. However, we are aware that it involves a lot of data, which is not always available, and leads to several computations due to the different preprocessing steps described in Section 3.2. For this reason, we tried to build an emulator with fewer inputs in X including only the variables from Table 3.2 at 700hpa and removing the solar and ozone forcings from Z. The results are reported in tables 3.3 and 3.4. They show that having more inputs increases the quality of the emulated series, but the “cheap” emulator presents satisfying results and can be considered as a serious option. For instance, in the first step of evaluation (Table 3.3), the cheap emulator provides a good representation of the spatial structure (S-Cor=0.94 vs 0.96 for the main emulator) but is more biased than the main emulator (-0.3 versus -0.18 °C). Having a specific selection of inputs for specific areas would prob-
ably increase the performance of the cheap emulator. However, this is not in the
spirit of the tool proposed here, which aims to be as simple and as general as possible.

Moreover, we also wanted to discuss the smoothing of the 2D inputs mentioned in
sections 3.2.2 and 3.2.5. There are two reasons behind it. First of all, it allows
getting closer to the effective resolution of the GCM models (Klaver et al, 2020).
Secondly, it deletes any local scale information that could remain after the upscaling
step when constructing the UPRCM inputs for the training (see Sect. 3.2.5). To
verify the usefulness of this step, we trained an emulator without prior smoothing
of the inputs. The results are reported under the name "No_Smoothing" in tables
3.3 and 3.4, in appendix C. The "No_smoothing" emulator performs slightly better
than the main emulator in perfect model world. It has for example a slightly better
RMSE (0.76°C versus 0.82°C for the main emulator) and a better ratio of variance
(96.2% versus 95.2%). However, in the GCM world, the "No_smoothing" emulator
shows less good results than the original emulator. The ratio of variance and the
spatial correlation over the different maps are the most adapted scores to compare
two emulators in GCM world (see sections 3.3.3 and 3.4.1). The ratio of variance of
the "No_smoothing" emulator falls to 83% for the worst points, while it is 95% for
the main emulator. The spatial correlation also confirms that it fails to reproduce
the entire spatial structure of the RCM. For example, when looking at the climate
change map, the main emulator's spatial correlation is equal to 98% while the one
of the No_smoothing emulator is 72%.

Finally, the last test on the inputs concerns the non-inclusion of the low resolution
surface temperature in the predictors. The motivation behind this choice comes
from the perfect model training detailed in Section 3.2.5. The intuition is that
too much high resolution information remains in TAS_UPRCM. The emulator uses
undoubtedly this information, leading to a less accurate downscaling of GCM fields.
To test this hypothesis, we trained an emulator including TAS, and the results are
referred as "TAS_included" in tables 3.3 and 3.4. The "TAS_included" emulator
performs better in the perfect model world but not in the GCM world. The most
striking score is probably the variance ratio. In GCM world, it ranges from 83%
to 120% which is closer to our very basic benchmark ([72; 128]%) than to the main
emulator ([95; 106]%). Moreover, the spatial correlations are also much less good
than for the main emulator. These results confirm clearly our hypothesis.

3.4.4 On the non reproducibility of neural network training

While neural networks have experienced considerable success over the last decades
and the number of applications is constantly increasing, they have also been largely
criticised for their lack of transparency due to an excessive number of parameters.
Several studies (see Guidotti et al, 2018, for review) have tried to provide the keys
to "open the black box". However, users of deep neural networks should also be
aware of their non-reproducibility. Indeed the training of deep neural networks with
GPUs involves several sources of randomness (initialisation, operation ordering,
etc.). A few recent studies have raised the issue for medical applications with
a UNet (Marrone et al, 2019; Bhojanapalli et al, 2021), but to the best of our
knowledge, no study has addressed it.

In order to document this issue for the RCM emulator developed in this paper and assess its robustness, we propose a Monte Carlo experiment where the same configuration of the emulator (as described in Section 3.2) is trained 31 times, resulting in 31 emulators. The results of this experiment are illustrated in Figure 3.14 and summarised in tables 3.3 and 3.4. The results of Section 3.3 are based on a randomly picked emulator, mentioned as the main emulator in appendix and plotted in darker green in Figure 3.14. All the emulators from the Monte Carlo experiment show consistent results with the ones presented in Section 3.3. The RMSEs (and correlated scores such as variance) are really close to each other ([$0.8; 0.86]°C$), which is expected since it is the loss function used to fit the neural network. However, the path taken to minimise the loss during training might vary from one emulator to another, leading to bigger differences in climatological metrics (see

<figure><img src="image_38.png" /><figcaption>**Figure 3.14:** Illustration of the results for the Monte Carlo experiment. (a) present a year time series for 4 given points and (b) the pdfs on the whole serie. The red line refers to the RCM truth, the dark green line is the main emulator, and the light green lines are the 30 emulators from the Monte Carlo experiment.</figcaption></figure>
tables). For instance, in the UPRCM evaluation step, the average error in the future
climatology varies from .003 to .139°C. The results from the GCM world evaluation
step (Table 3.4) are similar to those in the perfect model framework. This sta-
bility shows again the robustness of the emulator when using GCM large scale fields.

We believe that the readers of this study and any potential emulator users should be
aware of this characteristic of deep learning neural networks. However, it is worth
mentioning that this randomness in the training of the network does not impact the
key conclusions of the results (Section 3.3), which prove their robustness.

3.5 Conclusion

This study aims to explore a novel hybrid downscaling approach that emulates
the downscaling function of a RCM. That is to say learning the transformation
of the large-scale climate information into a local climate information performed
by a regional climate model. Here, we develop this approach for the near sur-
face temperature and a Southwest European domain. This new method, called
RCM-emulator, is designed to help increase the size of the high-resolution regional
simulation ensembles at a lower cost.

To achieve this overall goal, we develop a specific conceptual framework. The
emulator is trained using existing RCM simulations which allows it to learn the
large scale/local scale relationship in different climates and in particular in future
climate. Simply speaking, the general functioning of a RCM can be broken down
into a large scale transformation and a downscaling function. To focus on the
downscaling function, the emulator is trained in a perfect model framework,
where both predictors and predictand come from the same RCM simulation.
This framework implies to carefully prepare and select the predictors, to ensure
that no unwanted high resolution information remains in the training predictors.
The emulator takes daily large-scale and low-resolution information as input and
produces daily maps of the near surface temperature at the RCM resolution. It is
worth noting that the downscaling function likely depends on the RCM choice. So
the emulator developed in this study is RCM-dependent.

Technically speaking, the RCM-emulator is based on a fully convolutional neural
network algorithm, called UNet. A key point of the emulator is the substantial
computational gain regarding RCM computational cost. Training the emulator
used in this study took two hours on a GPU. While this time depends on the target
domain size it will never exceed some hours even for much bigger domains. Once
the emulator is trained, the downscaling of a new low resolution simulation takes
less than a minute. These time lengths do not include the preparation of the inputs;
nevertheless the gain remains evident when a RCM simulation involves weeks of
computation on a super-computer.

The emulator is evaluated in both perfect model and GCM worlds. The results
show that the emulator generally fulfils its mission by capturing very well the
transformation from low resolution information to the high resolution near surface
temperature. Firstly, the emulator is robust to different sources of input, which
validates our conceptual framework. Secondly, the emulator succeeds very well in reproducing the high resolution spatial structure and daily variability of the RCM. The perfect model evaluation shows that the emulator is able to reproduce the original series almost perfectly. Moreover, it appears clearly that training the emulator in future climate improves its ability to reproduce warmer climate. Nevertheless the emulator shows some limitations in accurately simulating extreme events and of the complete climate change magnitude. Future work should focus on these two aspects to further improve the emulator ability. It is worth noting that similar emulators could easily be built over different domains, for different RCMs and for surface variables other than near-surface temperature.

This UNet emulator has been compared to two other possible emulators based on more classical ESD methods using the well adapted perfect model framework. The exercise allows us to put into perspective the result of the UNet emulator and validate our choice as it outperforms these simpler emulators. Moreover, it is worth mentioning that this study does not target to identify the best possible emulator but aims to shed light on the feasibility to emulate the RCM complexity at high-frequency and high-resolution, and to propose one good technical solution to do it. We are fully aware that our method will certainly be one among others and we hope that the scientific community will propose new and probably better solutions for RCM emulators in the next years.

Finally, and even if it was not our original goal, this study highlights, as others before, the RCM-GCM inconsistencies at large scales. As the emulator focuses on the downscaling function of the RCM, it does not learn to reproduce these large scale transformations. This raises the question of how to evaluate the emulator when downscaling GCM simulation. An RCM run might not be the correct reference in such a case. Secondly, it also puts into question the final role of a RCM-emulator. Under the strong hypothesis that the large scale transformation carried out by the RCM results from physics or forcings inconsistencies – which will require further investigation – the emulator provides a high-resolution simulation that is corrected from the GCM-RCM large-scale inconsistencies.
## 3.6 Appendices

### 3.6.1 Emulators inter-comparison

In order to validate the neural network based emulator presented in this study, we proposed two other emulators relying on standard empirical statistical downscaling methods: CDFt and a Multiple Linear Regression. The idea is not to identify the best emulator, but to give a baseline and put into perspective the performance of the UNet architecture. The methods were chosen because both are well established and have been widely applied. Moreover they are good example of to two different families of ESD methods, Model Output Statistic (MOS, for CDFt) and Perfect Prognosis (PP, for MLR). Both methods are part of the large ESD method intercomparison project, VALUE (Gutiérrez et al, 2019).

CDFt was introduced by Michelangeli et al (2009) and used in several studies as Vrac et al (2012); Lavaysse et al (2012); Vigaud et al (2013); Famien et al (2018). It is based on the transformation of the Cumulative Distribution Function (CDF) of the low resolution variable into the CDF of the high resolution variable. In our case the low resolution variable is the smoothed surface temperature from the UPRCM conservatively interpolated on the RCM grid. We took the smooth TAS in order to be as coherent as possible with the UNet emulator inputs. The high resolution variable we try to reproduce is then the RCM near surface temperature TAS. To understand the theory behind CDFT, let $F_{T,Lr}$ be the CDF of the low resolution surface temperature from the training set, and $F_{T,Hr}$ the CDF of the high resolution surface temperature from the training set. Then CDFt relies on the assumption that there exists a transformation $T$ such that

$$T(F_{T,Lr}(x)) = F_{T,Hr}(x). \qquad (3.10)$$

Then by replacing $x$ by $F_{T,Lr}^{-1}(u)$, with $u$ any probability in $[0, 1]$ we get a simple definition of $T$:

$$T(u) = F_{T,Hr}(F_{T,Lr}^{-1}(u)). \qquad (3.11)$$

And under the assumption that $T$ is the same in the training set than in the evaluation set it comes that:

$$T(F_{E,Lr}(x)) = F_{E,Hr}(x), \qquad (3.12)$$

with $F_{E,Lr}$ (resp $F_{E,Hr}$) the CDF of the low (resp high) resolution surface temperature in the evaluation set. Then $F_{E,Hr}$ can be determined as:

$$F_{E,Hr}(x) = F_{T,Hr}(F_{T,Lr}^{-1}(F_{E,Lr})). \qquad (3.13)$$

It is then possible to reconstruct the high resolution series from the low resolution one, following a Quantile-Quantile approach between $F_{E,Hr}$ and $F_{E,Lr}$.

The second method is a Multiple Linear Regression (MLR), which has also been used in multiple studies (Huth, 2002; Huth et al, 2015; Manzanas et al, 2018; Joshi et al, 2015). It takes exactly the same list of inputs as the UNet emulator (Table 3.2), but by taking a single grid point for the 2D variables. The predicted high resolution temperature is then:
$$ \hat{Y} = \alpha_0 + \sum_{i=0}^{M} \alpha_i P_i, \qquad (3.14) $$

with M the number of predictors.

Both CDFt and MLR emulators are trained in the same conditions as the UNet

<figure><img src="image_39.png" /><figcaption>Figure 3.15: (a) Daily probability density functions from the RCM truth, the UNet, the CDFt and the MLR emulators at 4 particular grid points over the whole simulation period. (b) (resp. (c and d)) Maps of performance scores of the UNet emulator (resp. the CDFt and MLR) with respect to the RCM truth computed over the whole simulation period. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added.</figcaption></figure>
emulator. The simulation used for training is the ALADIN63 simulation forced
by CNRM-CM5, covering the 1950-2100 period with the RCP8.5 scenario from
2006. The ALADIN63 simulation forced by CNRM-CM5, covering the 2006-2100
period with the RCP4.5, is used to evaluate and compare the emulators. The
evaluation is made in perfect model framework as it gives the best conditions to
properly compare the results. The results are reported in Figures 3.15 and 3.16,
where the two simpler emulators are compared to the UNet Emulator used in the
study. Figure 3.15 is showing the performances of the three emulators regarding the
performance metrics defined in section 3.2.6. Both MLR and CDFt show a good
reproduction of the pdfs, the daily variance of the original RCM truth, as well as a
good temporal correlation. However regarding all these scores the UNet emulator
performs better, and this is well summarized by the RMSE maps, where the UNet
Emulator shows much more accurate results. Figure 3.16 presents the ability of
the emulators to reproduce some climatological statistics that are usually studied
in climate studies. The conclusions are very similar, the MLR and CDFt emulators
have good results. However the UNet emulator captures better the complexity of
the high resolution temperature simulated by the RCM. This results is clear when
watching at the spatial correlation of the Figure 3.16.

To conclude, both MLR and CDFt are reasonable candidates for the emulator, as
they were for statistical downscaling, however the UNet emulator provides better
results, specially in the more complex regions. Moreover, as a comparison, the very
simple multiple linear regression involves 233 thousand parameters over the whole
target domain when it uses a single grid point as input for each of the 2D variables
used. However, if we want to fit a MLR with the full 2D grid of input as for UNet,
the model would imply about 250 times more parameters (16*16), so more than the
UNet. For these reasons we decided to mainly present the UNet emulator.
<figure><img src="image_40.png" /><figcaption>**Figure 3.16:** (a) Maps of climatologies, (b) Number of days over 30°C and (c) the 99th percentile of daily temperature, for the present-climate period (2006-2025) and for the climate change signal (2080-2100 minus 2006-2025), for the RCM truth, the UNet, the MLR and CDFt emulators in perfect model mode. On each line, the two last maps show the error maps of the emulators. For each map, the spatial mean and super-quantiles (SQ05 and SQ95) are added, as well as the spatial correlation and spatial RMSE with respect to the RCM truth maps.</figcaption></figure>