La fonte de la neige qui s'accumule dans le Haut Atlas du Maroc représente une ressource en eau importante pour les habitants des vallées, car elle recharge les nappes souterraines et contribue au débit des oueds. Pourtant, l'état du manteau neigeux est mal connu à cause de l'absence de mesures in situ et de stations météorologiques à haute altitude. Pour répondre à ce défi, nous proposons dans cette thèse d'assimiler des données des surfaces enneigées dérivées des images à haute résolution spatiale et temporelle de la mission Sentinel-2 dans un modèle distribué de manteau neigeux forcé par des données de réanalyses climatiques. Cet objectif général nous a conduit à identifier dans un premier temps la résolution spatiale optimale pour simuler le manteau neigeux tout en limitant les temps de calcul. Ensuite, nous avons évalué l'apport des données de réanalyses météorologiques MERRA-2 pour estimer l'état du manteau neigeux dans un bassin versant peu étudié et important en terme de ressource en eau. Enfin, nous montrons l'implémentation et l'évaluation du schéma d'assimilation sur une saison complète dans un bassin versant pilote. Les résultats encourageants ouvrent la voie à une estimation régionale de l'équivalent en eau de la neige dans le Haut Atlas.

Introduction 1.1 Importance de la neige dans le Haut Atlas

La neige saisonnière fournit de l'eau à environ 2 milliards de personnes dans le monde [START_REF] Barnett | Potential impacts of a warming climate on water availability in snow-dominated regions[END_REF]. En particulier, dans les régions méditerranéennes et semi-arides, la fonte du manteau neigeux qui se forme chaque année dans les zones de montagnes est d'une grande importance hydrologique et économique [START_REF] Iglesias | Challenges to manage the risk of water scarcity and climate change in the Mediterranean[END_REF][START_REF] Viviroli | Mountains of the world, water towers for humanity : Typology, mapping, and global significance[END_REF][START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF]. Par ailleurs, dans ces mêmes régions, la pression sur la ressource en eau s'accentue. En effet, la demande en eau augmente à cause de la forte croissance démographique et du développement économique, ce qui accélère l'usage de l'eau dans tous les secteurs notamment en agriculture. C'est le cas dans la région de Marrakech au Maroc, où le développement socio-économique cause une demande drastique en eau (Agence Japonaise de coopération internationale, 2008). Alors que les précipitations annuelles moyennes dans les plaines sont faibles entre 150 mm et 200 mm, les précipitations atteignent 800 mm dans le Haut Atlas [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]. Une part importante de ces précipitations de haute altitude tombe sous forme de neige (entre 20% et 80% selon [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]). Ainsi, les montagnes accumulent en hiver un stock d'eau qui est précieux pour les régions en aval. Ce stock est d'autant plus important qu'il se déverse dans les oueds au printemps quand l'irrigation est indispensable pour soutenir la croissance des cultures. La contribution actuelle de la fonte des neiges au débit annuel des oueds atlasiques varie entre 15% et 45% selon (Boudhar et al., 2009) (Fig. 1.2).

Malgré l'importance du manteau neigeux comme ressource en eau au Maroc, la connaissance de ses propriétés physiques (e.g. hauteur de la neige, équivalent en eau, taux de fonte) reste pourtant très limitée. En effet, à cause d'un accès et d'une maintenance difficiles, les mesures in situ et les stations météorologiques d'altitude sont trop rares vis-à-vis de la forte variabilité spatio-temporelle du manteau neigeux en montagne [START_REF] Boudhar | Spatial distribution of the air temperature in mountainous areas using satellite thermal infrared data[END_REF]Baba et al., 2018a).

Dans le contexte marocain, peu de recherches se sont intéressées au suivi de l'enneigement. Une des premières études date de 1989 par [START_REF] Abidi | Martinec-Rango snowmelt runoff model applied to the Tillouguit Basin of Morocco[END_REF]. Elle consistait à modéliser les débits dans le bassin versant de Tillouguit (2554 km 2 ) en utilisant le modèle SRM (Snowmelt Runoff Model) [START_REF] Martinec | The snowmelt-runoff model (SRM) user's manual[END_REF]. Au début des années 2000, d'autres études sur l'enneigement dans les montagnes marocaines ont vu le jour. La première étude traitait de la sublimation et de la fonte de la neige dans la montagne de M'Goun à partir d'un jeu de données météorologiques et de données satellitaires [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF]. Les résultats de cette étude ont montré que la sublimation a un impact très important sur le bilan d'énergie du manteau neigeux dans les saisons froides et sèches, et contribue à 44% de pertes. Ensuite, dans le cadre du projet SUDMED [START_REF] Chehbouni | An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions : The SUDMED programme[END_REF] puis du LMI TREMA [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF], qui ont pour objectif de développer les connaissances nécessaires pour une gestion intégrée des ressources hydrologiques dans la région de Marrakech [START_REF] Chehbouni | An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions : The SUDMED programme[END_REF][START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF], plusieurs études se sont concentrées sur le versant nord du Haut Atlas dans le bassin de Tensift. Certaines études ont tiré partie des satellites d'observation de la Terre pour suivre l'évolution des surfaces enneigées [START_REF] Chaponniere | A combined high and low spatial resolution approach for mapping snow covered areas in the Atlas mountains[END_REF][START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Boudhar | Télédétection du manteau neigeux et modélisation de la contribution des eaux de fonte des neiges aux débits des Oueds du Haut Atlas de Marrakech[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF][START_REF] Ahmed | Outils spatiaux pour le suivide la dynamique de la neige sur l'ATLAS MAROCAIN : Changements climatiques et ressources en eau[END_REF]. D'autres études étaient consacrées à la modélisation de l'évolution du manteau neigeux à l'échelle de la station [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF]. Ces études ont permis de proposer une première évaluation de la contribution hydrologique du manteau neigeux et de caractériser la variabilité inter-annuelle des termes du bilan d'énergie du manteau neigeux : la contribution de la fonte de neige aux débits des rivières est importante et diffère d'un bassin à autre. En moyenne cette contribution dans les bassins du Tensift varie de 15% à 45%. De plus, il a été constaté que la durée d'enneigement entre février et mars avait tendance à augmenter à l'opposé de l'extension des surfaces enneigées qui avait tendance à décroître au cours de la période 2000-2013 [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF] ont aussi observé que la saison de fonte débute plus tôt qu'auparavant. Enfin, la pression sur la ressource en eau devrait s'aggraver nettement à l'avenir comptetenu de la réduction attendue du débit de printemps causée par le changement climatique [START_REF] Rochdane | Climate change impacts on water supply and demand in Rheraya watershed (Morocco), with potential adaptation strategies[END_REF][START_REF] Marchane | Climate change impacts on surface water resources in the Rheraya catchment (High Atlas, Morocco)[END_REF].

Néanmoins, avant cette thèse, il n'y avait pas d'étude permettant d'estimer la distri- bution spatio-temporelle de l'équivalent en eau du manteau neigeux dans le Haut Atlas. Pourtant c'est bien l'estimation du stock de neige dans les bassins versants qui intéresse les gestionnaires des ressources en eau comme l'agence du bassin hydraulique de Tensift, car cette information permet de mieux planifier la distribution de l'eau pour les irrigants en aval. Ce déficit de connaissance est lié au manque de données in situ d'une part et au fait que l'équivalent en eau (snow water equivalent en anglais, SWE) en montagne ne peut pas être estimé directement par satellite [START_REF] Dozier | Estimating the spatial distribution of snow water equivalent in the world's mountains[END_REF], un problème qui n'est pas spécifique aux montagnes marocaines.

Estimating the spatial distribution of snow water equivalent (SWE) has therefore been deemed the most important unsolved problem in snow hydrology.

Dozier et al. (2016)

Partant de ce constat, notre travail s'articule donc autour de cet objectif : comment estimer l'équivalent en eau de la neige dans le Haut Atlas ?

Pour y répondre, nous présentons dans la partie suivante les bases de la formation et l'évolution du manteau neigeux (Chap. 1.2) et des différentes approches pour son suivi : mesures in situ (Chap. 1.3.1), modélisation (Chap. 1.3.2) et télédétection (Chap. 1.3.3) afin de définir une stratégie générale pour répondre à cet objectif.

Évolution du manteau neigeux

Après chaque précipitation neigeuse, la neige s'accumule à la surface et subit plusieurs transformations au fil des jours en fonction de paramètres locaux qui sont liés à la topographie et aux effets météorologiques (e.g. température, précipitations, vent) à proximité de la surface. Ces effets atmosphériques sont résumés dans la figure 1.4. Pour évaluer l'évolution du manteau neigeux, le bilan de masse est utilisé. Ce bilan consiste à calculer le gain ou la perte de masse d'un manteau. Les précipitations neigeuses (P ) représentent un gain, tandis que la fonte (R), l'évaporation (E evap ) et la sublimation E sub sont des pertes. En l'absence du transport de la neige par le vent et par avalanches, l'équation de Figure 1.3 -Différentes formes de précipitations selon la température de l'air. Au niveau des hautes altitudes où la température est négative, les précipitations sont sous forme solide (cf. neige). Dans les zones où la température est élevée, les précipitations sont sous forme liquide. Alors que dans les zones où la température est proche de zéro, les deux différentes formes de précipitations sont possibles.

Chapitre 1 -Introduction l'évolution de la masse du manteau est exprimée comme suit :

dM dt = P -E sub -E evap -R (1.1)
Pour résoudre l'équation 1.1 il faut également résoudre l'équation du bilan d'énergie qui permet de déterminer l'énergie interne du manteau et ainsi la quantité d'énergie disponible pour la fonte. Elle s'écrit sous la forme suivante, en considérant que les flux vers l'extérieur du manteau sont comptés positivement :

dH dt = -(SW ↓ + SW ↑ + LW ↓ + LW ↑ + H s + H l + H p + G) (1.2)
-SW ↓ : Le rayonnement solaire, apporte de l'énergie principalement dans les longueurs d'ondes visibles et du proche infrarouge ([400 nm, 800 nm] -[800 nm, 2500 nm]). En absence d'absorption atmosphérique (ciel parfaitement clair) la valeur maximale sur une surface horizontale exposée perpendiculairement aux rayons du Soleil atteint 1361 W.m -2 . Ce flux est modulé par la valeur de l'angle solaire local qui dépend de la date et de la latitude, et par l'absorption atmosphérique. La quantité de rayonnement solaire reçue par le manteau neigeux dépend de la pente et de l'exposition. Le rayonnement solaire joue un rôle primordiale dans le bilan d'énergie du manteau neigeux [START_REF] Zuzel | Relative importance of meteorological variables in snowmelt[END_REF][START_REF] Marks | Climate and energy exchange at the snow surface in the alpine region of the sierra nevada : 1. meteorological measurements and monitoring[END_REF]. -SW ↑ : rayonnement solaire réfléchi par la surface du manteau neigeux. Si la neige est encore fraîche (fort albédo) le rayonnement sera réfléchi vers le ciel et donc peu absorbé par le manteau neigeux. Dans le cas inverse, le rayonnement absorbé jouera un rôle primordiale lors de la fonte [START_REF] Pomeroy | Physical properties of snow. Snow ecology : An interdisciplinary examination of snow-covered ecosystems[END_REF]. L'albédo de la neige est affecté par le dépôt d'impuretés, et par le métamorphisme qui diminue la surface spécifique de la neige [START_REF] Dumont | Optical remote sensing of snow cover[END_REF]. -LW ↓ : le rayonnement thermique infra-rouge [4 µm,100 µm] émis principalement par les couches d'air proches du sol. Il dépend de la température et de la colonne de vapeur d'eau. Si le ciel est nuageux, les nuages renvoient de la chaleur sous forme de rayonnement de grande longueur d'onde ce qui accélère la fonte. -LW ↑ : le manteau neigeux émet un rayonnement infra-rouge en fonction de sa température de surface. -H p : flux de chaleur originaire des précipitations. Il est important lorsque la pluie tombe sur le manteau neigeux, sinon il est négligeable. -Flux de chaleur sensible turbulent H s et flux de chaleur latente turbulent H l : correspondent respectivement aux flux de chaleur transportés par le mouvement de l'air et à la sublimation des cristaux de neige ou l'évaporation de l'eau liquide du manteau neigeux. Les équations de transfert simplifiées de ces flux sont dites bulk [START_REF] Male | Caractérisation de l'enneigement sur l'atlas marocain par le capteur MODIS et relation avec le climat (période 2000-2011)[END_REF]. Elles s'expriment comme suit :

H s = D H u z (T a -T s ) (1.
3)

H l = D e u z (e a -e s ) (1.4)
T est la température ( o )C. u z représente la vitesse du vent à 2 m. e représente la pression de vapeur. a et s dénotent l'air et la surface respectivement. D H et D E sont des coefficients dépendants de la stabilité et l'uniformité du terrain [START_REF] Pomeroy | Physical properties of snow. Snow ecology : An interdisciplinary examination of snow-covered ecosystems[END_REF]. Ainsi, le vent et la température de l'air et le gradient de l'humidité entre la surface du manteau neigeux et l'air au-dessus, jouent un rôle vital lors de la résolution du bilan d'énergie (Eq. 1.2). Par conséquent, ils contribuent fortement à la fonte. -Flux de chaleur vers le sol G : correspond aux échanges de chaleur avec le sol sous le manteau neigeux.

Chapitre 1 -Introduction

Effet de la topographie sur l'évolution du manteau neigeux

La topographie du terrain est un paramètre clé pour la détermination des variables de l'équation du bilan d'énergie (Eq. 1.2) comme résumé ci-dessous :

-L'altitude : Plus on monte en altitude, plus la pression diminue et par conséquent l'air se dilate ce qui le rend plus froid. Il est estimé que la température décroît selon un gradient qui est de l'ordre de -6.5 o C par kilomètre [START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF]. Ainsi les zones enneigées situées à basse altitude reçoivent un rayonnement atmosphérique et un flux turbulent de chaleur sensible plus important, ce qui accélère leur fonte. L'altitude joue aussi un rôle essentiel lors de la détermination des précipitations, puisque plus on monte en altitude, plus les précipitations sont plus importantes [START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF]. -Humidité relative : L'humidité relative est définie comme une mesure du rapport entre le contenu en vapeur d'eau de l'air et sa capacité maximale à en contenir. La quantité de vapeur d'eau que l'air peut contenir dépend de la température de l'air [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. Un air plus chaud peut contenir plus d'eau qu'un air froid. Ce qui explique qu'en général, l'humidité dans les basses altitudes est plus élevée que dans les hautes altitudes (où la température est trop basse (Chap. 1.2.2)). -Vent : Généralement la vitesse du vent augmente avec la hauteur dans la troposphère moyenne. Ainsi, le vent aux sommets est plus fort que le vent à basse altitude. Cette augmentation de vitesse du vent en fonction de la hauteur est en raison de la rapide diminution de la traînée du frottement avec la surface terrestre [START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF]. -L'orientation et la pente locale modulent l'intensité du rayonnement solaire reçu par le manteau. L'équation qui décrit ceci est :

SW ↓ = S(φ dir cos i + φ dif cos Z) (1.5)
S est la constante solaire, φ dir (W.m -2 ) représente transmissivité directe lorsque le ciel est dégagé, et φ dif (W.m -2 ) représente la transmissivité diffuse. Z est l'angle zénithale solaire, i est l'angle entre le rayon solaire direct et la surface inclinée, qui dépend de la pente locale et de son orientation par rapport à l'azimut solaire. -L'équation précédente ne tient pas compte des ombres projetées et du rayonnement diffus causés par le relief environnant [START_REF] Corripio | Vectorial algebra algorithms for calculating terrain parameters from DEMs and solar radiation modelling in mountainous terrain[END_REF][START_REF] Sirguey | Simple correction of multiple reflection effects in rugged terrain[END_REF].

Sublimation

Lorsque l'air est très sec et le ciel est dégagé (le rayonnement solaire incident est élevé), une partie du manteau neigeux "disparaît" par sublimation. Il faut 2838 kJ/kg à 0°C à pression standard pour que la neige se sublime.

Représentativité des forçages météorologiques

Comme mentionné plus haut, la température de l'air varie selon l'altitude. Mais en réalité, le gradient de température n'est pas fixe, et varie spatio-temporellement selon les conditions atmosphériques. [START_REF] Stahl | Comparison of approaches for spatial interpolation of daily air temperature in a large region with complex topography and highly variable station density[END_REF] ont comparé douze méthodes d'interpolation de température en Colombie Britannique au Canada. Ils ont remarqué que le manque des données au-dessus de 1500 m peut introduire des incertitudes importantes lors de la détermination du gradient de température. Dans le même sens, une étude préliminaire dans le bassin de la Rheraya au Haut Atlas marocain a montré que ce gradient variait entre 5.03 o C et 6.54 o C (Fig. 1.5).

Les précipitations ont tendance à augmenter avec l'altitude [START_REF] Barry | Mountain weather and climate[END_REF][START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF]. Mais à l'instar de la température, le facteur de correction de précipitation avec l'altitude varie en fonction du temps et de l'espace [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. Sa variabilité est même bien supérieure à celle de la température. Le gradient de précipitation peut parfois s'inverser avec l'altitude [START_REF] Scaff | Meteorological interpretation of orographic precipitation gradients along an Andes west slope basin at 30°s (Elqui Valley, Chile)[END_REF].

Dans tous les cas, les gradients météorologiques avec l'altitude sont des représentations très simplifiées de phénomènes climatiques complexes, et ils n'ont pour but que de décrire une variabilité spatiale au premier ordre.

Enfin, malgré son rôle important sur l'évolution et la redistribution du manteau neigeux, la distribution spatio-temporelle de la vitesse et de la direction du vent est encore mal connue en montagne [START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF]. De surcroît, les anémomètres dans les zones Chapitre 1 -Introduction montagneuses peuvent être confrontés à des vents violents ou peuvent être affectés par la neige et par la suite fourniront des résultats erronés.

Pour synthétiser, la détermination des forçages météorologiques dans les zones montagneuses est délicate à cause de leur forte variabilité spatio-temporelle. Par ailleurs, le nombre limité de stations météorologiques aux hautes altitudes pousse à extrapoler les mesures des régions localisées aux moyennes et basses altitudes, ce qui engendre de fortes incertitudes. C'est dans ce sens où [START_REF] Barry | Mountain weather and climate[END_REF] suggère que dans les zones montagneuses le réseau météorologique devrait être 10 à 20 fois dense que dans les zones de plaine afin d'assurer une estimation distribuée et précise de ces forçages. Malheureusement, l'application de ses recommandations est délicate (voir Chap. 1.3.1).

Suivi de l'évolution du manteau neigeux

La compréhension de l'évolution du manteau neigeux, sa description et ses caractéristiques impliquent l'usage de mesures in situ, de modèles de neige ou encore de données télé-détectées. Cette section présente ces différentes méthodes, leurs apports et limites.

Mesures in situ

Il existe plusieurs techniques in situ d'observations du manteau neigeux. Les capteurs ultrasoniques installés sur des station météorologiques automatiques déterminent la hauteur de la neige de façon continue. D'autre part, des campagnes de mesure de densité, de la profondeur et de l'équivalent en eau peuvent être effectuées le long d'un transect. Dans le Haut Atlas du Maroc, il n'existe qu'une seule station de mesure de la hauteur de neige pérenne installée à Oukaimeden (3230 m.a.s.l) 1 . Ce site de mesure a été sélectionné de telle façon qu'il soit accessible d'une part, et loin des pistes de ski et de randonnées d'autre part afin d'éviter l'influence humaine lors de la formation du manteau. Pour mesurer la hauteur de la neige, des relevés verticaux à l'aide d'un cylindre métallique sont effectués. A chaque point de mesure, trois carottes sont prélevées à ses alentours pour déterminer la profondeur moyenne du point initial. Afin de transformer ces hauteurs de neiges mesurées en quantité de l'équivalent en eau, nous déterminons sa densité (ρ mes ) en mesurant la masse de la calotte à l'aide d'une balance électronique. Une fois ces deux mesures sont disponibles, nous obtenons une estimation du SWE en utilisant la relation. 1.6 (Sturm et al., 2010).

SWE = HS ×

ρ mes ρ eau (1.6) 1. Il y a un capteur ultrasonique sur le toit du refuge du Neltner mais la représentativité de la mesure est discutable ρ eau est la densité de l'eau (1000 Kg.m -3 ). HS représente la profondeur de la neige (snow height).

Ces mesures aident à améliorer la compréhension de l'évolution du manteau neigeux à Oukaimeden, et donnent une vue globale sur l'enneigement du Haut Atlas, mais ne permettent pas de capturer l'hétérogénéité spatiale de la distribution de la neige et ne peuvent pas être spatialisées à l'échelle du bassin versant. Les campagnes aériennes LI-DAR ou la stéréoscopie par drone pourraient être une bonne solution pour résoudre la problématique de la variabilité spatiale du manteau neigeux. Pour l'instant le coût du LIDAR est prohibitif au Maroc, et les drones du commerce ont une couverture spatiale trop limitée pour être appliqués à l'échelle des bassins versants qui portent un intérêt pour la gestion de l'eau (>200 km 2 ).

Modélisation numérique de l'évolution du manteau neigeux Définition

La modélisation consiste à conceptualiser un système observé afin de le reproduire. Un modèle possède des paramètres fixes, des variables d'entrées qui sont transformées en sorties (variables simulées) (Fig. 1.7). En nivologie, la modélisation permet d'estimer les variables souhaitées en hydrologie (e.g. hauteur de la neige, équivalent en eau). Dans la littérature, il existe différentes approches pour cela. Pour simplifier, nous pouvons distinguer trois méthodes de simulation numérique de l'évolution du manteau neigeux :

-Approche ponctuelle (1D) : consiste à modéliser les caractéristiques du manteau neigeux en un point déterminé qui est généralement la position d'une station mé- Un point commun entre ces différentes approches est qu'elles sont très sensibles à la qualité des données météorologiques. En effet, une erreur dans les forçages se propagera au cours du temps et augmentera les incertitudes des sorties.

Télédétection du manteau neigeux

Définition

La télédétection peut être définie comme l'ensemble des techniques permettant d'acquérir des informations sur un objet à travers un instrument de mesure qui est loin de l'objet. En général, ce terme est lié à l'usage des satellites afin d'observer le globe terrestre. Cette acquisition consiste normalement à mesurer le rayonnement électromagnétique réfléchi ou émis dans un domaine de fréquence donné (e.g visible, proche infrarouge, micro-ondes).

Utilité

Dans les zones montagneuses l'accès est difficile notamment en hiver, ce qui complique la collecte de données sur le manteau neigeux, notamment dans les zones semi-arides méditerranéennes [START_REF] Boudhar | Spatial distribution of the air temperature in mountainous areas using satellite thermal infrared data[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. La télédétection optique (multispectrale) permet de remédier à ce déficit d'observations. La signature spectrale de la neige dans le domaine visible et moyen infrarouge (Chap. 1.3.3) permet de déterminer la surface enneigée. A l'heure actuelle, il existe plusieurs capteurs qui fournissent des images à différentes résolutions spatiales (de la résolution modérée 500 m -1000 m, aux hautes résolutions < 30 m) avec un temps de revisite acceptable (du journalier au 16 jours) pour le suivi du manteau neigeux (Dumont and Gascoin, 2017)(Tab. 1.3.3). Avec une revisite de 5 jours, une résolution de 20 m et une couverture globale, la mission Sentinel-2 a bouleversé la télédétection optique. Elle permet de suivre l'évolution de la surface enneigée des différentes régions du monde avec une grande finesse (Fig. 1.101.9.

Figure 1.9 -Illustration de l'apport de la télédétection spatiale pour le suivi de l'enneigement dans les zones montagneuses.

Contrairement à l'imagerie micro-onde, les produits de la télédétection optique sont affectés par les nuages. Dans le Haut Atlas marocain, un avantage est que la probabilité qu'un pixel soit couvert par un nuage est plus faible que dans les massifs sous climat tempéré comme les Alpes ou les Pyrénées. Néanmoins, la présence de nuages, combinée à la forte variabilité temporelle du manteau neigeux dans le Haut Atlas, incite à privilégier des missions ayant un temps de revisite court afin d'avoir un nombre important d'images sans nuages durant la saison de neige. C'est le cas des missions Sentinel-2 (5 jours), Venµs (2 jours), Aqua/Terra MODIS (1 jour).

Propriétés optiques de la neige et télédétection

La réflectance de la neige varie fortement avec la longueur d'onde. Dans le domaine visible, la neige est fortement réfléchissante. En revanche dans le moyen infrarouge, sa réflectance est très faible, ce qui permet de la distinguer des nuages qui ont généralement une réflectance élevée dans le visible et le moyen infrarouge (Fig. 1.11). Ainsi, l'identification du manteau neigeux est relativement aisée lorsque l'image satellitaire est multispectrale ou hyperspectrale et qu'elle contient des bandes dans le visible et le moyen infrarouge (MODIS, Landsat 8, Sentinel-2).

Extraction des surfaces enneigées à partir des images satellitaires

Il existe différentes méthodes de classification des surfaces enneigées à partir d'images satellitaires que nous résumons ci-dessous :

1.4. Problématique -Classification non-supervisée avec un indice spectral : L'indice de neige normalisé (Normalized Difference Snow Index en anglais, NDSI) se base sur le contraste entre la réflectance de la neige dans le visible et moyen infrarouge [START_REF] Dozier | Spectral signature of alpine snow cover from the Landsat Thematic Mapper[END_REF]. Le NDSI a été fréquemment utilisé pour déterminer les zones enneigées dans les images Landsat et MODIS [START_REF] Salomonson | Development of the aqua MODIS NDSI fractional snow cover algorithm and validation results[END_REF][START_REF] Hall | Normalized-difference snow index (NDSI)[END_REF] Il existe donc différents capteurs optiques permettant de suivre la surface enneigée jusqu'à une échelle très fine. Mais ces données ne permettent pas de suivre le taux de fonte, l'évolution de la hauteur de la neige ni celle de l'équivalent en eau, qui représentent les variables clés de l'hydrologie nivale. Des travaux récents ont montré qu'il est possible d'estimer la hauteur de la neige à partir d'un couple stéréoscopiques d'images à très haute résolution spatiale [START_REF] Marti | Mapping snow depth in open alpine terrain from stereo satellite imagery[END_REF]. Cependant l'application de cette méthode à l'échelle du Haut Atlas est loin d'être réalisable actuellement, compte-tenu du coût élevé des images stéréoscopiques (e.g 6.44E/km2 à travers AirBus.). D'autre part, les précisions atteignables (50 cm dans le meilleur des cas) sont relativement élevées par rapport à la hauteur de neige moyenne dans le Haut Atlas.

La télédétection passive dans le domaine des micro-ondes permet d'estimer la quantité de l'équivalent en eau à l'échelle de 25 km×25 km en inversant les températures de brillance [START_REF] Foster | Quantifying the uncertainty in passive microwave snow water equivalent observations[END_REF]. Mais la résolution spatiale des données micro-ondes s'avère être un facteur bloquant pour les zones montagneuses [START_REF] Dozier | Estimating the spatial distribution of snow water equivalent in the world's mountains[END_REF].

Problématique

Chacune des méthodes proposées auparavant (mesures in situ, modélisation et télédétection) présente des avantages et des limitations pour répondre à la problématique de

Synthèse de la problématique

Cette thèse vise à mettre en place une méthode qui permet de calculer l'état du manteau neigeux dans le Haut Atlas et en particulier le SWE sans donnée de terrain. Pour cela nous tâcherons de répondre aux questions suivantes :

-Quelle est la résolution spatiale permettant de simuler l'évolution du manteau neigeux en un temps raisonnable et sans dégrader l'information sur le SWE ? -Peut-on reconstituer l'enneigement dans le Haut-Atlas en utilisant des données de réanalyses météorologiques ? -Peut-on améliorer les résultats de ces simulations par assimilation de données satellitaires ? Dans les piémonts, la consommation de l'eau est en forte croissance à cause de la croissance démographique, du développement du tourisme (comme les terrains du golf) et surtout du développement de l'irrigation (Agence Japonaise de coopération internationale, 2008).

Une étude approfondie de [START_REF] Juillerat | Etude hydrologique comparative de 5 bassins versants du Haut Atlas marocain[END_REF] s'intéressait à l'analyse des données hydropluviométriques (Fig. 2.1) entre 1963 et 2003 à partir de données fournies par l'agence du bassin hydraulique du Tensift (ABHT). Nous reprenons ci-dessous les principaux résultats de cette étude.

Le climat est influencé par la proximité de l'océan Atlantique qui apporte un courant océanique froid [START_REF] Boudhar | Télédétection du manteau neigeux et modélisation de la contribution des eaux de fonte des neiges aux débits des Oueds du Haut Atlas de Marrakech[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. La moyenne des précipitations entre 1972 et 2002 dans le bassin versant de Tensift varie de 100 mm à 800 mm selon un gradient est-ouest lié à la topographie du bassin (Fig. Ainsi, la fonte de la neige dans les zones d'altitude du Haut Atlas contribue à l'alimentation de l'oued de Tensift et fournit de l'eau aux zones irriguées de la plaine du Haouz. [START_REF] Boudhar | Télédétection du manteau neigeux et modélisation de la contribution des eaux de fonte des neiges aux débits des Oueds du Haut Atlas de Marrakech[END_REF].

Bassins choisis

Dans cette thèse nous nous sommes focalisés sur les deux sous-bassins versants de l'Ourika et de la Rheraya (Fig. 2.1). Dans les deux cas, le régime hydrologique de ces bassins est fortement influencé par la fonte de la neige (Fig. 2.3). Le choix du bassin de la Rheraya est basé principalement sur la disponibilité des données météorologiques, et aussi parce que ce bassin a été l'objet de plusieurs études nivologiques lors des deux dernières décennies (par ex. [START_REF] Juillerat | Etude hydrologique comparative de 5 bassins versants du Haut Atlas marocain[END_REF][START_REF] Chaponniere | A combined high and low spatial resolution approach for mapping snow covered areas in the Atlas mountains[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF][START_REF] Richard | Evaluation d'un modèle pluie-debit sur le bassin versant de la Rheraya dans le Haut Atlas marocain[END_REF]. Le bassin de l'Ourika est en revanche caractérisé par un déficit de données et d'études nivologiques, malgré un apport hydrologique bien plus significatif (Fig. 2.3). Il représente donc un bassin de choix pour appliquer une méthode de modélisation du manteau neigeux qui n'utilise pas des données in situ.

Ces deux sous-bassins sont situés au centre Ouest du Maroc (Fig. 5.4). Leurs caractéristiques physiographiques sont présentées dans le tableau. 2.1.1. Le Toubkal est situé au sud du bassin de la Rheraya. Les deux bassins de notre étude reçoivent un taux de précipitations relativement important pour une région semi-aride (Fig. 2.2). Cela est lié à leurs altitudes élevées (Tab. 2.1.1). Les débits de ces bassins sont marqués par une forte variabilité inter-annuelle [START_REF] Juillerat | Etude hydrologique comparative de 5 bassins versants du Haut Atlas marocain[END_REF]. différentes caractéristiques de ces stations hydro-météorologiques.

Données disponibles

Modèle utilisé

SnowModel a été développé principalement par Glen E. Liston à l'Université du Colorado (Liston and Elder, 2006b,a;[START_REF] Liston | Simulating complex snow distributions in windy environments using SnowTran-3D[END_REF]. C'est un modèle d'évolution du manteau neigeux à base physique qui résout les équations du bilan de masse et d'énergie de façon distribuée (sur une grille régulière). Cette représentation permet de prendre en compte l'effet de la topographie et de la végétation dans une modélisation spatialisée. SnowModel est alimenté par des forçages météorologiques mesurés par des stations météorologiques (précipitation, température et humidité relative de l'air, vitesse et direction du vent à 2 m, rayonnement solaire incident de courte longueur d'onde et rayonnement atmosphérique de grande longueur d'onde) à un pas de temps horaire à journalier. Il prend aussi comme entrée le modèle numérique de terrain (MNT), la carte de l'occupation de sol, ainsi qu'un fichier de configuration qui définit les options et les paramètres du modèle (voir Chap.3.2.4,Tab. 3.3). SnowModel permet de simuler plusieurs variables décrivant l'évolution du manteau neigeux (e.g. hauteur de la neige, équivalent en eau, fonte). Snow-Model a été utilisé dans plusieurs études nivologiques en zones montagneuses [START_REF] Hiemstra | Snow redistribution by wind and interactions with vegetation at upper treeline in the Medicine Bow Mountains, Wyoming, USA[END_REF][START_REF] Hiemstra | Observing, modelling, and validating snow redistribution by wind in a wyoming upper treeline landscape[END_REF][START_REF] Mernild | Climatic control on river discharge simulations, Zackenberg River drainage basin, northeast Greenland[END_REF][START_REF] Gascoin | Wind effects on snow cover in Pascua-Lama, Dry Andes of Chile[END_REF][START_REF] Szczypta | Impact of climate and land cover changes on snow cover in a small Pyrenean catchment[END_REF].

Ce modèle de manteau neigeux est composé des quatre sous-modèles ci-dessous (Fig. 2.7).

-MicroMet : modèle de distribution des variables météorologiques dans la zone de l'étude. Il se base sur la méthode de Barnes [START_REF] Barnes | A technique for maximizing details in numerical weather map analysis[END_REF][START_REF] Koch | An interactive barnes objective map analysis scheme for use with satellite and conventional data[END_REF] [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF] :

T 0 = T stn -τ (Z 0 -Z stn ) (2.3)
T stn ( o C) est la température observée à l'altitude Z stn (m). T 0 ( o C) est la température à l'altitude de référence Z 0 (m) (Z 0 = 0 m). τ ( o C)m -1 est le facteur de correction de température (temperature lapse rate en anglais). Ensuite l'interpolation de Barnes est appliquée. Enfin, le modèle numérique terrain (MNT) et l'équation 2.4 sont utilisés pour ajuster le niveau de référence de température aux altitudes données par le MNT :

T (x, y) = T 0 -τ (Z -Z 0 ) (2.4)
T (x, y) est la température à l'altitude Z du point de grille (x, y) donnée par le MNT.

La spatialisation des précipitations est légèrement différente par rapport à celle de la température. En premier temps, les précipitations sont interpolées à la grille du modèle en utilisant le schéma d'analyse de Barnes. Ensuite, une surface de référence est générée à partir des altitudes des stations. Dans un dernier temps, les précipitations modélisées P (x, y) (mm) au point (x, y) à l'altitude Z (m) sont calculées comme suit :

P (x, y) = P 0 [ 1 + χ(Z -Z 0 ) 1 -χ(Z -Z 0 ) ] (2.

5)

P 0 représente les précipitations interpolées à la première étape à l'altitude interpolée Z 0 . τ (mm.m -1 est le facteur de correction de précipitations (precipitation lapse rate en anglais). Un exemple de spatialisation des précipitations mettant en exergue l'effet de l'altitude sur sa distribution est présenté dans la figure 2.9. L'estimation du rayonnement solaire est effectuée à travers un ensemble d'équations qui prennent en compte l'influence de la couverture nuageuse σ c , les transmissivités directe φ dir et diffuse φ dif ainsi que la pente et l'orientation du pixel. Premièrement, le coefficient de fraction nuageuses (σ c ) est déterminé en utilisant la valeur de l'humidité relative à 700 hpa (RH 700 ) [START_REF] Walcek | Cloud cover and its relationship to relative humidity during a springtime midlatitude cyclone[END_REF][START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF] :

σ c = 0.832 × exp[ RH 700 -100 41.6 ] (2.6)
Ensuite, les transmissivités directe φ dir et diffuse φ dif sont déterminées à partir des équations 2.7 et 2.8 [START_REF] Burridge | The Meteorological Office operational 10-level numerical weather prediction model[END_REF] : où T wb représente la température humide (wet bulb temperature en anglais) qui est à son tour estimée à partir de la température de l'air [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]). -SnowTran-3D : ce module permet de calculer la redistribution du manteau neigeux en utilisant l'effet du vent [START_REF] Liston | Simulating complex snow distributions in windy environments using SnowTran-3D[END_REF]. Nous ne détaillerons pas son fonctionnement, puisque nous avons désactivé ce module soit par manque de données de direction du vent, soit parce que son temps de calcul est trop élevé.

φ dir = (0.6 -0.2 × cosZ)(1 -σ c ) (2.7)

Préparation des données 2.3.1 Modèle numérique terrain

Une connaissance détaillée de la topographie est importante pour le calcul des forçages. Pour cela, nous avons généré un MNT de la Rheraya à 4 m de résolution à partir de trois couples stéréos panchromatiques (6 images) issus du capteur du satellite Pléiades.

La chaîne de traitement des couples stéréo-Pléiades se base principalement sur des commandes fournies à l'outil Ames Stereo Pipeline (ASP), un outil développé par la NASA, conçu pour le traitement des images stéréoscopiques captées par des satellites. Les étapes de cette chaîne sont présentées dans le chapitre. 3.2.3. La figure. 2.10 synthétise ces processus. [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF].

Cartes de surfaces enneigées

Saison 2015-2016

En Juin 2015, Sentinel-2A a été mis en orbite. Depuis cette date, il fournit des images tous les 10 jours. Après le lancement du Sentinel-2B en Mars 2017 la résolution temporelle de Sentinel-2 est rapidement passée à 5 jours. Les images de ce satellite constituent une mine d'or puisqu'elles ont à la fois une bonne résolution spatiale et temporelle, ce qui permet de bien suivre les surfaces enneigées. Pour cette thèse, nous nous sommes limités à l'usage des données de Sentinel-2A pour la saison 2015-2016 à cause de l'absence de données de validation hydro-météorologiques dans notre zone d'étude après cette année.

Autres saisons

Dans le chapitre 4 nous avons mené une étude sur l'évolution du manteau neigeux du bassin d'Ourika depuis 2000 jusqu'au 2018. Les images à haute résolution spatiale couvrant toute cette période n'ont pas une résolution temporelle adéquate pour suivre l'évolution du manteau neigeux (seules les données Landsat-7 et Landsat-8 avec un temps de revisite de 16 jours sont disponibles). Nous avons donc utilisé les produits neige du satellite MODIS (avec un temps de repassage journalier). La chaîne de traitement de ces données pour interpoler les pixels masqués par les nuages est décrite dans [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF].

Outils d'évaluation

Afin d'évaluer les simulations effectuées dans les différents chapitres, plusieurs métriques sont utilisées.

Évaluation des débits

A l'exutoire des deux bassins versants (Rheraya et Ourika) l'Agence de Bassin Hydraulique du Tensift (ABHT) relève quotidiennement les hauteurs d'eau. Des jaugeages sont idéalement réalisés au pas de temps mensuel à cause des contraintes logistiques. Le jaugeage se base sur la technique de l'exploration des champs de vitesse [START_REF] Boudhar | Télédétection du manteau neigeux et modélisation de la contribution des eaux de fonte des neiges aux débits des Oueds du Haut Atlas de Marrakech[END_REF]. Cette fréquence de jaugeage affectera la qualité des débits observés notamment après des pluies violentes.

La simulation des débits à partir des sorties du modèle n'est pas directe, puisque que le modèle de neige ne calcule pas l'évapotranspiration ou l'infiltration, seulement la lame d'eau de ruissellement (runoff) issue des pluies et de la fonte de la neige. Donc, pour définir la qualité des simulations, nous nous sommes basés sur la comparaison entre les débits et le runoff du modèle à un pas de temps assez large et de façon qualitative seulement.

Évaluation des hauteurs de la neige

Une station météorologique est installée à Oukaimeden (Fig. 2.12), à la limite des deux bassins versants Rheraya et Ourika. Elle comporte une sonde ultrasonique qui fournit en continu des mesures de la hauteur de la neige. Ces données sont moyennées au pas journalier et utilisées pour qualifier les simulations. Notons que pour quelques années les mesures de la hauteur de la neige sont erronées ou manquantes. De ce fait, les années d'observations de la hauteur de la neige se sont limitées sur celles où les données sont d'une bonne qualité (2004-2005, 2005-2006, 2008-2009, 2010-2011, 2013-2014 et 2014-2015).

Évaluation des surfaces enneigées Conversion SWE-SCA

Dans le but d'évaluer les sorties du modèle avec des données de surfaces enneigées par satellite, une conversion SWE-SCA (snow water equivalent-snow cover area) doit être effectuée. Dans le cadre de ce travail, nous nous sommes concentrés seulement sur deux types de SDC. Une description de chacune d'elles est présentée ci-dessous :

-Utilisation d'un seuil sur l'équivalent en eau, à partir duquel on considère le pixel comme complètement enneigé (SCA=1) et en dessous, complètement déneigé (SCA=0) (Fig.

Matrice de confusion

De nombreuses études utilisent l'évolution de surfaces enneigées intégrées sur tout un bassin pour valider un modèle de neige. Pour tirer partie de l'information spatialisée dans les produits de télédétection, une méthode couramment utilisée en télédétection consiste à calculer la matrice de confusion entre les deux images que l'on compare (Fig. 2

.14).

Dans notre étude, la matrice de confusion est composée de 4 cellules (Fig. 2.14). Les valeurs de ces cellules sont déterminées à partir d'une comparaison pixel par pixel de la grille des surfaces enneigées modélisées avec celles observées par télédétection.

-TP (true positive en anglais) : le modèle prédit que le pixel est enneigé et c'est vrai d'après les observations. -TN (true negative en anglais) : le modèle prédit que le pixel est non-enneigé et c'est vrai d'après les observations. -FP (false positive en anglais) : le modèle prédit que le pixel est enneigé alors que c'est faux d'après les observations. -FN (false negative en anglais) : le modèle prédit que le pixel est non-enneigé alors que c'est faux d'après les observations.

Des métriques sont associées à la matrice de confusion permettant de synthétiser les performances du modèle.

Précision

Cet indice permet de déterminer la capacité du modèle à reproduire correctement les différentes classes (classe de neige et classe sans neige). Il est donc définit comme la somme des bonnes estimations (TP et TN) sur l'ensemble de toutes les estimations (TP, TN, FP (2.12)

Plus le HSS est proche de 1 plus l'estimation est bonne. Un HSS proche de 0 correspond aux mauvaises simulations du SCA. La formule Eq. 2.12 est valable seulement lorsqu'aucune classe n'est prédominante par rapport à l'autre [START_REF] Notarnicola | Snow cover maps from MODIS images at 250 m resolution, part 2 : Validation[END_REF] Baba, M. W., Gascoin, S., Kinnard, C., Marchane, A., and Hanich, L. Effect of digital elevation model resolution on the simulation of the snow cover evolution in the High Atlas (preprint https: // rdcu. be/ bIJGa

Motivations et résumé de l'étude

En théorie, la distribution spatiale du SWE dans les zones montagneuses est fortement corrélée à la topographie. En effet, les précipitations, la température, le rayonnement solaire et le vent sont dépendants de l'altitude, de l'orientation et des pentes du terrain (Sec. 1.2.2). L'usage d'un MNT à très haute résolution spatiale devrait donc permettre de mieux représenter la distribution des forçages météorologiques, et par suite l'évolution du manteau neigeux. D'un autre côté, le temps d'exécution d'un modèle de neige distribué augmente rapidement avec la résolution du MNT. Ainsi, avant de réaliser une simulation sur une longue période temporelle ou bien de mettre en place un schéma d'assimilation qui implique la génération d'un grand nombre de simulations, il était crucial de répondre à cette question : Dans le contexte du Haut Atlas, quelle est la résolution spatiale du MNT permettant de préserver l'effet de la topographie tout en s'exécutant dans un temps opérationnel ? Afin de répondre à cette problématique, nous avons préparé un MNT à 4 m de résolution à partir de trois couples stéréo Pléiades (0.5 m de résolution spatiale). Ce MNT a été rééchantillonné à 8 m, 30 m, 90 m, 250 m, 500 m et 1000 m. Ces différents MNT ont été utilisés en entrée du modèle de neige SnowModel afin de simuler l'évolution du manteau neigeux dans le bassin de la Rheraya à différentes résolutions. La qualité des résultats des différentes simulations (8 m à 1000 m) a été évaluée à travers deux jeux de données : des mesures in situ de la hauteur de la neige à Oukaimeden (3230 m.a.s.l.) et des surfaces enneigées dérivées d'une série temporelle d'images à haute résolution du satellite Formosat-2 (8 m de résolution) acquise spécifiquement sur ce bassin au cours de la saison 2008-2009 dans le cadre du programme SUDMED [START_REF] Chehbouni | An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions : The SUDMED programme[END_REF]. Cette série était encore unique en son genre avant l'avènement de la mission Sentinel-2 en 2017 avec le lancement de Sentinel-2B. L'objectif de cette étude est d'évaluer la performance du modèle distribué du manteau neigeux SnowModel aux différentes résolutions des MNTs (de 8 m à 1000 m). Les résultats obtenus indiquent que le modèle a une bonne performance de 8 m à 250 m et la qualité diminue à partir de 500 m. L'analyse des sorties du modèle suggère que cette dégradation est due au lissage de la topographie. En plus de l'évaluation des simulations, une analyse complémentaire du MNT de la zone d'étude a montré que le semivariogramme de la pente se rompt à partir de 200 m. Ce constat est en harmonie avec les résultats de simulations. Nous avons conclu que la résolution de 250 m était suffisante pour assimiler les SCA dans un modèle de neige pour estimer le SWE dans cette zone.

L'étude a été menée dans le bassin pilote de la Rheraya, mais nous spéculons que les résultats peuvent être transposés à d'autres bassins du Haut Atlas puisque les caractéristiques climatiques et physiographiques seront similaires. Cependant, cela devra être fait avec précaution, puisque notre étude se concentre seulement à la saison 2008-2009. De surcroît, le module SnowTran-3D n'a pas été activé à cause de la non-disponibilité des mesures de vent.

Article : Effect of digital elevation model resolution on the simulation of the snow cover evolution in the High Atlas

Abstract

The snow melt from the High Atlas represents a crucial water resource for crop irrigation in the semi-arid regions of Morocco. Recent studies have used assimilation of snow cover area (SCA) data from high resolution optical sensors to compute the snow water equivalent (SWE) and snow melt in other mountain regions. These techniques however require large model ensembles and therefore a challenge is to determine the adequate model resolution, which yields accurate results with reasonable computation time. Here we study the sensitivity of an energy-balance model to the resolution of the model grid for a pilot catchment in the High Atlas. We used a time series of 8 m resolution SCA maps with an average revisit time of 7.5 days to evaluate the model results. The DEM was generated from Pléiades stereo-images and resampled from 8 m to 30 m, 90 m, 250 m, 500 m and 1000 m. The results indicate that the model performs well from 8 m to 250 m but the agreement with observations drops at 500 m. This is because significant features of the topography were too smoothed out to properly characterize the spatial variability of meteorological forcing, including solar radiation. We conclude that a resolution of 250 m might be sufficient in this area. This result is consistent with the shape of the semivariogram of the topographic slope, suggesting that this semivariogram analysis could be used to transpose our conclusion to other study regions.

Key Points :

-A distributed energy-balance snow model is applied in the High Atlas for the first time.

-The model performance decreases at resolution coarser than 250 m. -This result is consistent with the semivariogram of the topographic slope.

Introduction

In many semi-arid and Mediterranean catchments, snow melt runoff is an essential water resource for farmers who need to irrigate their crops during the dry summer season [START_REF] Viviroli | Mountains of the world, water towers for humanity : Typology, mapping, and global significance[END_REF][START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF]. In the Tensift region of Morocco near Marrakech, the major rivers have their headwaters in the High Atlas mountains. The contribution of snow melt to their mean annual discharge was estimated to range from 15 % to 45 % [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]. However, this contribution is probably more important from a water resource management perspective because (i) the annual flood due to snow melt occurs at this period of the year when irrigation is crucial to the crop growth, (ii) the melt water flood is more stable and easier to use than the flash floods triggered by heavy rainfalls, which often cause more harm than good in this region [START_REF] Zkhiri | Regional flood frequency analysis in the High Atlas mountainous catchments of Morocco[END_REF].

The analysis of optical remote sensing data showed that the snow cover in the High-Atlas is characterized by high spatial and temporal variability [START_REF] Chaponniere | A combined high and low spatial resolution approach for mapping snow covered areas in the Atlas mountains[END_REF][START_REF] Boudhar | Long-term analysis of snow-covered area in the Moroccan High-Atlas through remote sensing[END_REF][START_REF] Male | Caractérisation de l'enneigement sur l'atlas marocain par le capteur MODIS et relation avec le climat (période 2000-2011)[END_REF]. However, these studies did not address the variability of the snow water equivalent (SWE) and snow melt, which are the main variables of interest for hydrologists and water managers. The variability of the SWE and snow melt was analyzed at the point scale using energy balance models forced by meteorological station data [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF][START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF][START_REF] López-Moreno | Different sensitivities of snowpacks to warming in Mediterranean climate mountain areas[END_REF]. The important results of these studies are that the main source of energy for melting is net radiation and that snow sublimation can be significant but with a high inter-annual variability. Computed sublimation ranged between 7 % and 25 % of the total ablation at Oukaimden (3230 m) on the windward slope of the High Atlas, while it reached 45 % at Tichki (3260 m) on the drier, lee-side of the High Atlas. These point-scale modeling studies, however, did not allow the computation of melt rates at the catchment scale. To our best knowledge, the only attempts to compute the snow melt at catchment scale in the High Atlas were based on degree-day models [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Boudhar | Apport des données FORMOSAT2 à la modélisation du contenu en eau du manteau neigeux du Haut Atlas marocain[END_REF]). Yet, the aforementioned point-scale modeling studies suggest that distributed energy-balance models should be more suitable to capture the processes driving the evolution of the snowpack in the High Atlas area since the standard degree day approach does not explicitly include the radiation and neglects sublimation. A drawback of distributed energy-balance models is that they require more input forcing data than simpler, empirical models. In the mountainous areas of Morocco, the in situ network of meteorological stations is extremely sparse, while the spatial variability of the climate is exacerbated by the complex topography of the terrain. The errors caused by the larges uncertainties in the input data can spread over the course of a snow season due to the temporal persistence of the snowpack properties and cause significant biases in the SWE and snow melt (Durand et al., 2008b).

The recent surge of high-resolution remote sensing data, e.g. the distribution of the entire Landsat archive and the advent of the Sentinel-2 mission, have opened new avenues to the field of snow hydrology in semi-arid and Mediterranean regions [START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF]. These data are free and enable to map the snow cover extent with a spatial resolution (20 m to 30 m) which is compatible with the variability of the snow cover in semi-arid regions. For instance, it enables to capture the large differences in snow cover area due to the contrasted melt rates between north and south facing slopes (Pimentel et al., 2017a). A limitation to the Landsat mission for snow cover applications is the low frequency of observation (Landsat-8 mission has a revisit time of 16 days). However, with Sentinel-2 the revisit time will be 5 days globally. In addition, the cloud cover in semi arid regions is more limited than in temperate regions. Thus the main limitation is rather that Landsat-8 and Sentinel-2 are optical sensors and hence do not allow retrieving the internal snowpack properties such as the SWE and melt rates. This is why these data must be combined with a snowpack model through the calibration of the model parameters or within a data assimilation framework which accounts for the model and observation errors. This approach is considered superior to approaches based on passive microwave remote sensing to estimate the SWE in semi-arid and Mediterranean mountain regions [START_REF] Dozier | Estimating the spatial distribution of snow water equivalent in the world's mountains[END_REF]. For instance, the assimilation of MODIS and Landsat SCA enabled to reduce the bias in the simulated SWE by 49 % in the Rio Grande headwaters (USA) (Durand et al., 2008b). More recently, the assimilation of Landsat SCA into a distributed energy balance model was used to reconstruct the SWE in the California Sierra Nevada and the subtropical Andes [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF]. The updated SWE estimation reduced the bias between modeled SWE and snow course by 60 % to 82 %.

These studies suggest that the assimilation of high resolution SCA maps in a distributed snow cover energy-balance model is a promising approach to improve our knowledge of the snow resource in the High Atlas. However, we lack perspectives on the applicability of this method since a snow cover energy-balance model has never been applied in the Atlas mountains, even without data assimilation (i.e. in an open loop). In particular, an important question before implementing a distributed data assimilation scheme is what is the adequate model spatial resolution ?

In theory, the highest model resolution should provide the best results since the spatial variability of the SWE in mountains is strongly linked to the topographic variability [START_REF] Molotch | Estimating the spatial distribution of snow water equivalent in an alpine basin using binary regression tree models : the impact of digital elevation data and independent variable selection[END_REF]. Key factors are the dependence of precipitation and temperature to elevation, the effect of slope and aspect on incoming solar radiation and the deflection of the wind field by the terrain [START_REF] Barry | Mountain weather and climate[END_REF][START_REF] Jost | The influence of forest and topography on snow accumulation and melt at the watershed-scale[END_REF]. In addition, higher resolution DEM may more accurately reflect late-lying snow patches than lower resolution DEM simulations. Melt water from these snow patches can be critical to sustain streamflow during the summer [START_REF] Freudiger | Snow redistribution for the hydrological modeling of alpine catchments[END_REF].. Previous numerical experiments suggest that the mean SWE simulated by a distributed snowpack can decrease by 10% by changing the model grid resolution from 25 m to 1000 m [START_REF] Schlögl | Sensitivity of Alpine3D modeled snow cover to modifications in dem resolution, station coverage and meteorological input quantities[END_REF]. In the High Atlas, our field experience suggests that the model should at least resolve the hillslope scale, i.e. a spatial scale with an order of magnitude of 10 2 m (Fig. 3.1).

While higher spatial resolution are expected to better capture the variability of the snow cover, this is at the expense of the computational cost. Evaluating the computation time of the model is crucial in the perspective of data assimilation, which requires a large number of model runs. [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF] used a model grid size of 90 m based on literature values [START_REF] Blöschl | Scaling issues in snow hydrology[END_REF][START_REF] Winstral | Long-term snow distribution observations in a mountain catchment : Assessing variability, time stability, and the representativeness of an index site[END_REF], but did not assess the impact of this choice on their results. Recently, [START_REF] Baldo | Implementation of a physiographic complexity-based multiresolution snow modeling scheme[END_REF]) developed a multi-resolution approach to optimize the computational cost, while preserving the necessary resolution for fine-scale snow processes. This scheme is based on a physiographic complexity metric snow modeling scheme and was evaluated with respect to a reference simulation at high-resolution. However, the results of the new scheme were not compared to observations. Other strategies (than changing the model grid resolution) were proposed to account for terrain heterogeneity on topographic-driven meteorological forcing while keeping low computational requirements. These approaches can be classified in two categories : (i) the subgrid approach [START_REF] Müller | A grid-and subgrid-scale radiation parameterization of topographic effects for mesoscale weather forecast models[END_REF][START_REF] Essery | Scaling and parametrization of clear-sky solar radiation over complex topography[END_REF][START_REF] Gagnon | A gibbs sampling disaggregation model for orographic precipitation[END_REF] ; (ii) the semi-distributed approach [START_REF] Younas | A strategy to represent impacts of subgrid-scale topography on snow evolution in the Canadian land surface scheme[END_REF][START_REF] Revuelto | Meteorological and snow distribution data in the Izas experimental catchment (Spanish Pyrenees) from 2011 to 2017[END_REF]. The first approach was used to represent forcing variables as probabilistic distributions (instead of a single mean value) at the grid cell scale. This approach was mainly developed for land surface schemes [START_REF] Arola | Effects of subgrid spatial heterogeneity on GCMscale land surface energy and moisture fluxes[END_REF][START_REF] Essery | Scaling and parametrization of clear-sky solar radiation over complex topography[END_REF], motivated by the fact that these models are designed to provide lower boundary conditions to regional and global climate models, which have kilometer resolution at best. To our best knowledge, however, there is no study which put together all forcing variables (radiation, wind, temperature, etc.) as subgrid processes to run a snowpack energy balance model in alpine terrain. The second option aims to explicitly resolve the effect of topography on pre-defined classes of topography. It has been notably used by Météo-France for operational avalanche forecasting and was recently implemented in the Canadian Land Surface Scheme [START_REF] Younas | A strategy to represent impacts of subgrid-scale topography on snow evolution in the Canadian land surface scheme[END_REF]. However, this approach is not practical to assimilate re-mote sensing observations which are distributed on a regular grid. In that perspective, it is more straightforward to run the model on a distributed grid. This also permits to explicitly represent snow redistribution processes such as the wind transport or avalanches [START_REF] Revuelto | Meteorological and snow distribution data in the Izas experimental catchment (Spanish Pyrenees) from 2011 to 2017[END_REF].

The objective of this study is to determine an adequate model resolution to simulate the snow cover evolution in the semi-arid context of the Atlas mountains. We focus on a snow-dominated pilot catchment in the Atlas mountains, where forcing and validation data are available. We adopted a practical approach, i.e. we sought the lowest resolution which did not deteriorate the model performances. According to [START_REF] Blöschl | Scaling issues in snow hydrology[END_REF] the processes that influence the snowpack variability across the landscape have correlation lengths ranging from 1 m (wind drift) to 10 km (climatic conditions), while the effect of solar radiation at hillslopes of different aspects has a typical process scale of 100 m. Hence, we generated digital elevation models (DEMs) with resolutions of 8 m, 30 m, 90 m, 250 m, 500 m and 1000 m of the catchment using Pléiades stereoscopic imagery. These DEMs were used to assess from which resolution the model deteriorates too much with respect to validation data. The simulations were done with SnowModel (Liston and Elder, 2006a), a distributed snow evolution model that was already tested in various regions, including semi-arid mountains (e.g. [START_REF] Gascoin | Wind effects on snow cover in Pascua-Lama, Dry Andes of Chile[END_REF]. Validation data include continuous snow depth measurements from an automatic weather station, but the model evaluation mainly relies on a time series of 8 m resolution SCA maps derived from Formosat-2 acquisitions during the 2008-2009 snow season. These Formosat-2 data can be seen as a surrogate of the Sentinel-2 data (the operational revisit frequency of 5 days was achieved in the Atlas during the 2017/2018 snow season). The accuracy of the simulations is balanced with the computation time. Finally, we investigate the consistency of the model results with the spatial variability of the topography in the study area, which allows us to propose a potential method to help determine the adequate model resolution to run an energybalance model beyond this particular location. To our best knowledge this study is the first to study the effect on model grid resolution based on in situ and high resolution remote sensing observations in a semi-arid mountainous environment. It is also the first evaluation of a distributed energy balance model in the High Atlas.

Study area and Data

Study area

We focus on the Rheraya catchment located in the High Atlas range (31.18 o N, 7.86 o W, Fig. 3.1). The catchment area is 220 km 2 and its elevation ranges from 1000 m to 4167 m asl (Toubkal, highest peak in North Africa). This catchment is characterized by a semi-arid climate, with irregular precipitation. Snow accumulation generally occurs above 1500 m asl between December and March and the melt season ends in June (Boudhar et al., 2007a) . However, melt events are not uncommon in winter [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. The Rheraya catchment was selected because it is a pilot site of the Tensift observatory, which is coordinated by the joint international laboratory TREMA [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF] within the framework of French-Moroccan cooperation. A number of meteorological stations were set up in the Rheraya catchment since 2003, including two high elevation stations (see Section 3.2.3). These stations provide critical data in this region where the meteorological network is otherwise extremely sparse (Sect. 3.2.2). In addition, the Rheraya river is by itself an important tributary of the Tensift river, which runs through the even more arid Haouz plain and the city of Marrakech (currently about 1M inhabitants in the metropolitan area, Fig. 3.2). However, only a very low fraction of the runoff reaches Marrakech since most of the runoff from the headwater catchments such as the Rheraya feeds the Haouz aquifer by infiltration in the river bed or is diverted for irrigation between the mountains foothills and the city. The Haouz aquifer is the main source of water for irrigation and human consumption in the region.

Data

We used hourly meteorological data from four Automatic Weather Stations (AWS). These stations provide standard meteorological variables (Tab. 3.1), and the snow height was recorded at Oukaimeden AWS using an acoustic snow gauge [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF].

All these stations were installed during the early phase of the SUDMED project between 2002 and 2004 [START_REF] Chehbouni | An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions : The SUDMED programme[END_REF] and were maintained as part of the Tensift observatory by the staff of the joint international laboratory TREMA [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF].

Only Oukaimeden station is equipped with a wind-sheltered Geonor T-200B precipitation gauge. The other stations (Aremd, Imskerbour, Tamatert and Asni) have tipping bucket rain gauges. However, these stations rarely receive snow in contrast to Oukaimeden station. The instruments at Oukaimeden were described and by [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF] who used the data to evaluate the energy fluxes at the station. The station is located near the summit of a peak and just outside of a ski run on a gentle north-facing slope (supplement Fig. S1).

A DEM of 4 m resolution was produced for this study specifically. It was derived from Pléiades stereoscopic imagery. The Pléiades system is a constellation of two identical satellites (named Pléiades 1A and 1B). Both platforms are on the same near-polar sun- synchronous orbit at 694 km of altitude and can acquire images at 70 cm resolution in the panchromatic band (480-830 nm). Due to its agility Pléiades can acquire two images of the same area with different angles along the same orbit. These stereoscopic images can be combined to derive a DEM. Pléiades spatial resolution is initially equal to 70 cm, but the images are provided at 50 cm after pre-processing.

We obtained three Pléiades stereo-pairs covering the Rheraya catchment. These stereopairs were acquired on 04-08-201504-08- , 17-08-201504-08- and 18-08-2015 (Fig. 3.3) (Fig. 3.3). The acquisition dates were constrained to summer to avoid the presence of snow. The images of 18-08-2015 are partly covered by clouds near their western edge. The base to height ratios (B/H) were respectively 0.32, 0.35 and 0.3 for each date. B/H ratios lower than 0.4 are recommended for mountainous terrain [START_REF] Hasegawa | DEM accuracy and the base to height (b/h) ratio of stereo images[END_REF]. The operator (Airbus Space and Defense) provides a Rational Polynomial Coefficient (RPC) file with each image to allow the geo-localization of the images. These images and associated RPC files were used to generate a high resolution DEM with the Ames Stereo Pipeline (ASP, version : 2. 5), an open source stereogrammetry software developed by NASA [START_REF] Shean | An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery[END_REF], which supports Pléiades data and is efficient in mountainous regions [START_REF] Marti | Mapping snow depth in open alpine terrain from stereo satellite imagery[END_REF]. We applied the following workflow based on ASP command line utilities :

-Alignment : in mountainous region classical correlation which consists to match analogue pixels based on the RPC georeferencing only may fail due to sharp relief. Hence we pre-projected the right and left Pléiades images on the SRTM1 DEM before aligning each pair separately using the map_project command. -Stereoscopy : this step consisted on generating point clouds from stereoscopic pairs using the stereo command.

-Rasterization : the three point clouds were rasterized to grids of 4 m resolution through the point2dem command. -Mosaicing and blending : the three DEMs were merged using the dem_mosaic command. -Gap-filling of cloud covered area : In cloudy areas we obtained no-data values because the stereo correlation failed in the previous step. We filled these gaps by merging the Pléiades DEM with the SRTM DEM. This was done by resampling the SRTM DEM to 4 m of resolution by cubic interpolation. Then, we merged both DEMs using the dem_mosaic command again but with a priority blending length option which gives priority to the Pleiades DEM. -Orthoimage : Once the DEM is complete we used it to orthorectify the panchromatic images at 1 m resolution in the WGS-84 UTM 29N spatial reference system. This orthoimages will thereafter allow us to validate horizontal accuracy of the DEM. -Filtering and gap-filling of artifacts : using a hill-shaded version of the DEM and a by loading the DEM in a 3D visualization application, we identified obvious artifacts (areas of noisy terrain and spurious spikes) in the DEM in the north-west of zone 1 near the cloudy area and in the south of zone 3 near the Toubkal peak. We manually delineated these erroneous pixels in a GIS software and changed their value to No-data (They represent 2 % of the data). These remaining gaps were isolated and covered areas generally smaller than the SRTM pixel size, hence we filled them using the gdal_fillnodata.py utility with a radius of 11 pixels as a maximum search distance (44 m). The final DEM is shown in Fig. 3.2.

To evaluate the DEM accuracy, 7 ground control points (GCP) were collected with a SX Blue dGPS in the Rheraya catchment on 14 May 2015. The mean absolute error in elevation at the GCP is equal to 4.72 m and the Normalized Median Absolute Deviation (indicator of the random error) is 4.10 m (defined as N M AD = 1.4826median(|∆h j -m ∆h j |) where ∆h j refers to the individual errors and m ∆h denotes the median of errors [START_REF] Höhle | Accuracy assessment of digital elevation models by means of robust statistical methods[END_REF]). Although higher vertical accuracies can be attained with Pléiades [START_REF] Berthier | Glacier topography and elevation changes derived from Pléiades sub-meter stereo images[END_REF][START_REF] Marti | Evolution of Ossoue Glacier (French Pyrenees) since the end of the Little Ice Age[END_REF] these results are within the specifications given by Airbus Defense and Space (vertical accuracy inferior to 10 m). To further assess the horizontal accuracy, we used GPS points collected during a 40 km long road trip with the Sx-Blue device mounted on the roof of the car. The Sx-Blue was set up with an acquisition frequency of 1 point per second (tracking mode). We superposed these points to the Pléiades ortho-image to verify that the points are aligned with the road (Fig. 3.4). We also digitized several roads subsections in different regions from the ortho-image. We found that 90 % of the control points collected along these road stretches are within the digitized road polygons. We estimate that roads in this area have a width of 3 to 8 m [START_REF] Boulaajoul | Modelling the snow distribution at two high arctic sites at Svalbard, Norway, and at an alpine site in central Norway[END_REF]. The snow cover area maps were derived from a time series of Formosat-2 images and have a spatial resolution of 8 m. Formosat-2 was launched in 2014 for the purpose of monitoring natural resources, forestry and environment [START_REF] Chern | On-orbit soh monitor of a satellite-a case study of FORMOSAT-2[END_REF]. Its spatial resolution is respectively 2 m for the panchromatic band and 8 m for the multispectral bands which were used in this research to extract the SCA. 19 Formosat-2 images were acquired from February 2008 to June 2009 with a constant view angle over the snow-dominated part of the Rheraya catchment (see the images extent in Fig. 3.5). 15 of them are cloudfree and exploitable to derive SCA, which corresponds to an average satellite overpass of 7.5 days, spanning a period starting with an almost completely snow covered scene until complete snow ablation. Images were geometrically corrected by cross-correlation with a SPOT reference image, which was previously geo-referenced in the UTM-29N projection system [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. Atmospheric correction was performed with the method of [START_REF] Hagolle | Correction of aerosol effects on multi-temporal images acquired with constant viewing angles : Application to Formosat-2 images[END_REF] for the retrieval of aerosol optical thickness. Then, the snow-covered masks were computed from a supervised classification based on bands 2 (Green), 3 (Red) and 4 (Near-infrared) using the parallelepiped approach as implemented in the ENVI software [START_REF] Lillesand | Remote sensing and image interpretation." minimum distance to means classifier[END_REF]. The training data were sampled in every image separately using three classes : snow, no snow and shadowed snow (between five and ten polygons of each class for each image). The two snow classes were then merged in one snow class. The images were visually inspected and if necessary additional training samples were collected to improve the classification. About 10 samples were used for each date and each class. The accuracy of the results was estimated to exceed 90% for all dates. The only modification made to the data presented in [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF] is that we refined the geolocalization of the Formosat-2 images stack using the Pléiades ortho-image (Sect. 3.2.3). Indeed we found a small offset between both datasets in the upper part of the catchment. The correction was done by identifying control points that were visible in both the Pléiades ortho-image and true color composites of the snow-free Formosat-2 images. These reference points were used to fit the parameters of a translation vector, which was then applied to the full Formosat-2 stack. Given the accuracy of the Pléiades ortho-image geolocalization that was estimated above (Sect. 3.2.3), we estimate that the Formosat-2 snow maps of this study have an absolute registration accuracy lower than 8 m. More importantly, the relative co-registration with the Pléiades DEM is lower than 7 m. This time series of high-resolution SCA maps is unique in the High Atlas. However, it is a precursor of the type of data that will be routinely available when the Sentinel-2 mission will be fully operational (5 days revisit time, cloud permitting). As such, it represents an excellent test bed to test and develop methodologies that will rely on Sentinel-2 data.

Land cover

The catchment land cover was obtained by a supervised classification of Google Earth Imagery and field observations [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. The surface was classified in three classes : bare soil (68 %), shrubs (29 %) and deciduous forest (3 %). The land cover map is used in SnowModel to determine various parameterizations linked to the vegetation effects on the snowpack. Here the vegetation has a very marginal effect on the simulations since most of the study domain is classified as bare soil.

Methodology

Model description and setup

We used SnowModel, a distributed snowpack evolution model (Liston and Elder, 2006a). SnowModel is composed of four sub-models : MicroMet, EnBal, SnowPack and SnowTran-3D. MicroMet generates a distributed meteorological forcing over the model domain from AWS measurements [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. The SnowPack submodel uses the outputs from MicroMet and the energy-balance submodel EnBal to simulate the evolution of the height and density of the snowpack (Liston and Elder, 2006a). We used the one-layer option as recommended by the default settings. The SnowTran-3D submodel redistributes the snow depth over the terrain depending on wind direction and speed (Liston et al., 2007a). However, we did not activate this module for this study due to the lack of in situ wind direction measurements.

The core of the model is based on the resolution of the surface energy balance of the snowpack as implemented in EnBal (Liston and Elder, 2006a).

Q m = (1 -α)Q si + Q li + Q le + Q h + Q e + Q c (3.1) Where Q si (W.m -2 ) is the incoming shortwave radiation, Q li (W.m -2 ) is the incoming longwave radiation, Q le (W.m -2 ) is the emitted longwave radiation, Q h (W.m -2 ) is the turbulent exchange of sensible heat, Q e (W.m -2 ) is the turbulent exchange of latent heat, Q c (W.m -2 ) is the conductive energy flux, Q m (W.m -2
) is the the energy flux available for melt and alpha is the surface albedo. We modified some aspects of MicroMet. Spatial interpolation in MicroMet is based on the Barnes objective analysis scheme, which uses a Gaussian distance-dependent weighting function to interpolate the station data to the regular model grid. Prior to this interpolation, elevation-dependent corrections are applied to remove the effect of elevation on the precipitation, air temperature and humidity, i.e. as if the domain was flat. Once the Barnes interpolation is performed, the data are lapsed back to the grid elevation as given by the DEM. In the case of air temperature, the elevation correction is a linear function of elevation :

T x = T stn -τ (Z x -Z stn ) (3.2)
Where T stn ( o C) is the observed station air temperature at Z stn elevation. T x ( o C) refers to the temperature in a given elevation at Z x .τ is a monthly lapse rate ( o C.km -1 ) given in Tab. 3.2.

To distribute the precipitation over the domain MicroMet interpolates the observed precipitation by using Barnes objective scheme. Then, it generates a topographic reference surface by interpolating also the elevation of stations. This topographic reference is used to adjust the precipitation with a non-linear function of elevation difference ( Eq. 4.3 ). The modeled precipitation P x (mm), at grid point x of elevation Z x , is equal to :

P x = P 0 × [1 + χ(Z x -Z 0 )] [1 -χ(Z x -Z 0 )] (3.3)
Where P 0 (mm) is the interpolated precipitation, Z 0 (m) is the interpolated station elevation. χ (km -1 ) is the precipitation correction factor, also given at the monthly timestep (Tab. 3.2). [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF] provides default values for air temperature lapse rates and the precipitation correction factors. However, these values were computed for our study area based on AWS observations. The air temperature lapse rate was computed by ordinary linear regression using the monthly averages (Table 3.2). The monthly precipitation correction factors were computed by least-square optimization of (Eq. 4.2).

MicroMet simulates the incoming shortwave radiation to the surface (Q si ) by taking into account the effect of the geometry of the local slope and assuming that the effects of the surrounding topography are negligible (e.g. no casted shadows).

Q si = S(φ dir cos i + φ dif cos Z) (3.4)
where S is solar beam radiation (1370 W.m -2 ) [START_REF] Kyle | The status of the Nimbus-7 earth-radiation-budget data set[END_REF], Φ dir (W.m -2 ) is the direct solar radiation, Φ dif (W.m -2 ) is the diffused solar radiation, Z is the solar zenith angle, i is the angle between direct solar radiation and a sloping surface.

cos i = cos B cos Z + (sin B sin Z cos(µ -ξ s )) (3.5)
B is the terrain slope, µ is the solar azimuth and ξ s the terrain slope azimuth (with south having zero azimuth). For more details the reader can refer to [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF].

When available, MicroMet assimilates shortwave measurements using a direct insertion method : the difference between simulated and observed shortwaves at the meteorological stations is calculated and interpolated to the model grid using the Barnes interpolation scheme ; the difference grid is then added to the simulated shortwave radiation. If a meteorological station is located on a grid cell with a high slope value, the simulated incoming shortwave at this station will be affected by the cos B factor. However, incoming shortwave radiation is typically measured in the horizontal plane by pyranometers. For this study it can be problematic since the terrain of the Rheraya catchment is very rugged. For example, Oukaimeden AWS is located on a ski hill. In addition, the local slope can vary significantly depending on the DEM resolution, which could introduce a spurious dependence of the radiation forcing to the input DEM resolution. Therefore, we modified the routine for radiation data assimilation as follows : first, we compute the incoming shortwave radiation for a flat domain. Then, we compute the differences between the simulated and measured shortwave radiation at the stations. These differences are distributed using the Barnes analysis and added to the shortwave radiation grid which is computed using the DEM. This first order correction assumes that the atmospheric clearness does not depend on the topography.

The wind speed, which influences the turbulent heat exchange, is also modified according to the topography (Liston et al., 2007a), but we did not change the associated parameters given that only two AWS measured the wind speed and wind direction was not recorded during the study period.

SnowModel was ran at an hourly time step from 01 November 2008 to 10 June 2009. All other parameters were set to default values (Tab. 3.3). We performed five simulations by changing only the model resolution (8 m, 30 m, 90 m, 250 m, 500 m and 1000 m). The Pléiades DEM was resampled using bicubic interpolation.

The resampling resolutions were chosen based on the following considerations : 8 m is the native resolution of Formosat-2 ; 30 m is the highest resolution that can be theoretically reached from standard DEMs (SRTM, ASTER) ; 90 m is the highest resolution of current snow reanalyses [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF] ; 250 m is the highest resolution at which the snow cover area can be extracted from MODIS (Notarnicola et al., 2013a) ; 500 m is the resolution of the standard MODIS snow product.

Evaluation of the simulations

The simulations were evaluated by comparing daily snow depth averages measured at the Oukaimeden AWS, with the simulated snow depth extracted at the corresponding grid cell. The root mean square error (RMSE) was used as model verification metric.

We also spatially evaluated the simulations using the Formosat-2 SCA maps. The validation domain was defined as the intersubsection of the catchment area and the extent of the Formosat-2 acquisitions. We used the simulated SWE to compute the snow presence or absence for each pixel in order to compare the simulations with the binary Formosat-2 snow maps. This conversion of SWE to SCA was done using a constant threshold, i.e. for each pixel, if SWE > SWE 0 , then SCA = 1, otherwise SCA = 0 (Fig. 3.6). We did not use a more complex snow depletion curve formulation [START_REF] Pimentel | Subgrid parameterization of snow distribution at a Mediterranean site using terrestrial photography[END_REF][START_REF] Kolberg | Interannual stability of grid cell snow depletion curves as estimated from MODIS images[END_REF][START_REF] Magand | Introducing hysteresis in snow depletion curves to improve the water budget of a land surface model in an Alpine catchment[END_REF] because these parameterizations introduce more parameters that we are not able to determine in our study area. In addition, parameters of the snow depletion curve are expected to change with the model resolution. To account for the uncertainty in the value of the SWE threshold we generated the simulated snow maps from a list of 50 equally-spaced SWE 0 values from 1 mm to 120 mm w.e. This range was chosen based on previous comparison between remote sensing snow maps and station data [START_REF] Klein | Validation of daily MODIS snow cover maps of the Upper Rio Grande River Basin for the 2000-2001 snow year[END_REF][START_REF] Gascoin | A snow cover climatology for the Pyrenees from MODIS snow products[END_REF]. All simulated snow maps were then resampled to the Formosat-2 resolution (8 m) using the nearest-neighbour method (Fig. 3.6). This enabled to compute a confusion matrix between each simulated snow map and the corresponding Formosat-2 snow map for each of the 15 acquisition dates.

From the confusion matrix we computed the Heidle skill score (HSS), a statistical index that is recommended when a class is largely dominant [START_REF] Notarnicola | Snow cover maps from MODIS images at 250 m resolution, part 2 : Validation[END_REF] (our domain is dominated by the snow-free class) : HSS = 2(T P × TN -FP × FN)/[(T P + FP) × (FP + TN) + (T P + FN) × (FN + TN)], where TP (true positive) is the number of cells which were classified as snow in the simulation and the Formosat-2 maps, FP the number of false positive, etc. When a class is dominating, it can be shown that : lim d→∞ HSS = 2TP/(2TP + FP + FN) (Notarnicola et al., 2013a). The perfect simulation has an HSS equal to 1 while the worst has an HSS equal to 0.

Model calibration

It is important to reduce as much as possible the bias in the model inputs to be able to pinpoint the effect of the model grid resolution. For instance, [START_REF] Baldo | Implementation of a physiographic complexity-based multiresolution snow modeling scheme[END_REF]) used a meteorological dataset that was optimized for snowpack modelling from Landsat observations. In our study area, the main source of bias presumably comes from the uncertainty in the distributed precipitation field. We used precipitation data from five stations for a catchment area of 220 km 2 . It is by far the highest density of precipitation gauge in the Atlas mountains. However, there remain large uncertainties on the high elevation precipitation since the elevation of the highest gauge is 3230 m asl, while the maximum elevation in the catchment is 4167 m asl (Tab. 3.1). Therefore, we have used the time series of Formosat-2 images to optimize the monthly precipitation correction factors in MicroMet (7 parameters from November to May included). The best correction factors were selected in the interval (0 -0.35 km -1 ) by Monte-Carlo sampling using the Heidle skill score over the entire study period as optimization variable. The correction factors interval was chosen based on [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. The model was initialized with the correction factors derived from the weather stations (Sect. 3.2). Given the large number of model runs that is required for such optimization, we have performed the optimization The simulated SWE at every Formosat-2 acquisition date was converted to binary snow maps using a SWE 0 threshold. Then, the simulated snow maps were resampled to the same grid as Formosat-2. This enables to compute a confusion matrix between both images, from which we can derive the statistics described in Sect. 3.2.4.

Results

Snow cover area Fig. 3.7 shows the evolution of the snow cover area in the Rheraya catchment before and after the optimization of the monthly precipitation lapse rates. It shows that the optimization has been successful to reduce a general bias in the model since the evolution of the snow cover area is closer to the Formosat-2 observations. Then, the simulated SCA maps were compared on a pixel basis to the Formosat-2 observations after resampling all the simulated snow maps to 8 m (Sect. 3.2.4,Fig. 3.6). This procedure was repeated for every Formosat-2 date and every SWE 0 . Fig. 3.8 shows the statistical distribution of the resulting HSS values at each resolution. The differences between the models are low from 8 m to 250 m, but the median of the HSS values decreases at 500 m and even more at 1000 m.

We selected the SWE 0 threshold which gives the highest HSS for each resolution and we computed the average of these SWE 0 values. The resulting threshold (40 mm) was used to compute the SCA on every date of Formosat-2 acquisition (15 dates) for all the resolutions. Then, we summed for each pixel the number of days with snow. The same calculation was made with the Formosat-2 data. This comparison shows how coarsening the model resolution modifies the snow patterns at the catchment scale. Simulations from 8 m to 250 m produce similar snow cover duration patterns. At 500 m and 1000 m the resolution is too coarse to capture the variability in the snow cover duration that reflects important terrain features like the main valley bottoms. The impact of the resolution is much less visible when the snow cover area is considered as an aggregated value at the catchment scale (Fig. 3.7), however the root mean squared error (RMSE) in the SCA is the largest at the coarsest resolutions 500 m and 1000 m (Tab 3.4).

Snow depth

The comparison with the in situ snow depth record shows that the model captures main dynamics of the snowpack at Oukaimeden, although there is a consistent overestimation of the ablation rates. The snow depth evolution suggests that this excessive ablation rate was partly compensated by an overestimation of the precipitation. However, simulations from 8 m to 250 m are closer to the observations than the 500 m and 1000 m simulations, which exhibit a large negative bias (Fig. 3.10 and Tab. 3.4). The large bias at 500 m and 1000 m can be attributed to the smoothing effect of DEM resampling on the elevation of high peaks. Indeed, Oukaimeden AWS is located on an isolated peak (supplement Fig. S1), therefore the DEM coarsening tends to reduce the elevation at this specific location. The elevation of the Oukaimeden grid cell ranges between 3238 m and 3231 m in the DEMs at 8 m to 250 m resolution, while it is only 3132 m and 2848 m at 500 m and 1000 m respectively.

Discussion

Changing the model resolution essentially means that we test the model sensitivity to the topographic resolution. Other model parameters may have a much larger impact on the model performance, however the model resolution has a special interest since it is the parameter that has the strongest impact on the computation time and memory usage. To better understand the effect of the model resolution on the simulations, we have plotted The effect of the model resolution on the energy balance is more evident in on the spatial distribution of the snow cover. We focused on the southern part of catchment, a high elevation region where the topography is highly variable with contrasted slopes and aspect. In this area we expect that the meteorological forcing are properly represented given the presence of the Neltner station (Tab. 3.1). As an example, we selected the 01 April to map the SWE at each resolution in this area, and we compared it to the Formosat-2 snow cover area (Fig. 3.12). It shows that the model is able to capture the spatial variability of the snow cover at the hillslope scale. The SWE is lower in southfacing slopes than in north-facing slopes. This spatial variability is not inherited from the precipitation forcing since the precipitation variability at this scale in the model is driven by the elevation (Eq. 4.3). Fig. 3.13 shows that the simulated snow patterns in this area are similar to the incoming shortwave radiation. At 250 m and below, the model is able to capture the fact that the slopes of the east-west oriented valleys are not equally snow-covered. Above 250 m the model does not capture anymore this slope effect. This observation can be verified on other acquisition dates during the melt season. This can be explained by the fact that the DEM resampling modifies the slope and aspect of the model grid cells, which are used in the calculation of the radiation fluxes (Eq. 3.4). As a Based on these observations, we computed the empirical semivariograms of the DEM slope in the study area. Here, we used the ASTER global DEM V2 rather than the Pléiades DEM, because it is the highest resolution DEM which is available globally at no cost and thus it could be used to do the same semivariogram analysis in other regions. We find that the semivariogram clearly flattens out at length greater than 500 m (Fig. These values are consistent with our previous findings regarding the effect of the resolution on the snow cover simulation. Indeed, the semivariogram analysis implies that above a grid spacing of approximately 200 m, the variability of the slope is not captured anymore by the DEM. Because the variance of the slope is the main parameter controlling the mean and variance of solar radiation over topographic surfaces [START_REF] Dubayah | Topographic distribution of clear-sky radiation over the Konza Prairie, Kansas[END_REF][START_REF] Mckenney | Calibration and sensitivity analysis of a spatially-distributed solar radiation model[END_REF][START_REF] Pohl | Spatial-temporal variability in solar radiation during spring snowmelt[END_REF][START_REF] Essery | Scaling and parametrization of clear-sky solar radiation over complex topography[END_REF], it could also explains why the model performance degrades significantly above 250 m in our case.

Conclusion

The objective of this study was to evaluate the performance of a distributed snowpack energy balance model at different resolutions. We applied SnowModel at different resolutions using AWS measurements as input data. We found that the model performed relatively well at resolutions from 8 m to 250 m but the performance dropped significantly at 500 m and 1000 m. The same conclusion was drawn from both in situ observations (snow depth at Oukaimeden station) and remote sensing data (Formosat-2 snow maps).

Since the Rheraya catchment can be considered as a representative catchment in the snow dominated areas of the High Atlas, these results suggest that a resolution of 250 m might be sufficient to simulate the SWE at the scale of the whole High Atlas range, especially in the perspective of a data assimilation scheme based on a particle filter. The reduction of computation time is substantial from 8 m to 250 m (Tab. 3.5). Current data assimilation scheme typically involve ensembles of 100 simulations [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF]. Below 250 m, we estimate that the computation times are actually prohibitive given our computer resources for an ensemble-based data assimilation experiment of several snow season at the larger scale of the High Atlas mountains.

The differences in the model outputs are the results of multiple interactions between the DEM and the equations used in MicroMet to spatially interpolate the precipitation, air temperature, humidity, wind speed, longwave and shortwave radiation. Thus, the results of this study are dependent on this model. However, the main equations used in MicroMet are commonly used to spatially interpolate the meteorological forcing in other models, hence it is possible that our conclusions remain valid with a different model. We interpret that the main difference are due to :

-The aggregation of the DEM to a coarser resolution tends to reduce the elevation of the highest peaks and to increase the elevation of valley bottoms. The net impact of this DEM smoothing on the snow accumulation is not obvious, however, since it will depend on the mean elevation of the rain-snow transition. -The aggregation of the DEM to a coarser resolution tends to decrease the slope of the model grid cells and can change the slope orientation (aspect). This modifies the radiative forcing and therefore the snow cover patterns during the melt season.

The results of this study must be taken with caution since there remain large uncertainties in the model input and output, and also because some topography-dependent processes were not taken into account such as the transport of the snow by the wind (due to a lack of wind observations) or avalanches (due to a lack of model parameterization). However, preliminary tests, using wind data from a climate model reanalysis, indicated that the simulated wind transport is low in the study area. Indeed, in SnowTran-3D the wind transport occurs when the surface friction velocity exceeds a threshold value, which depends on the snow density [START_REF] Liston | A snow-transport model for complex terrain[END_REF]. In the High Atlas, the snowpack Table 3.5 -Computational time of the simulations at each spatial resolution when the model is run on a 800 Mhz processor with 16 Go of memory.

Spatial resolution 8 m 30 m 90 m 250 m 500 m 1000 m Computing time (1 run) 15 h 1 h 8 mn 1 mn 12 s 3 s temperature stays close to 0 o C, melt events are common during the accumulation season and thus, the surface snow density increases rapidly, which inhibits the wind transport.

Another limitation of the model is the lack of a casted shadow algorithm. We estimate that the additional fraction of the catchment area that would be in the shadow if such algorithm was implemented is 23% at 10 :00, 4% at 12 :00 LT, 0.1% at 14 :00, 0.3% at 16 :00 and 21% at 20 :00 local time. Further work is necessary to evaluate the impact of this deficiency on the snow cover simulation. More advanced radiation schemes are available to improve the spatial interpolation of the radiation in complex terrain [START_REF] Aguilar | Topographic effects on solar radiation distribution in mountainous watersheds and their influence on reference evapotranspiration estimates at watershed scale[END_REF]. Finally, an important limitation is that the study is limited to a single snow season, while the study area is characterized by a strong inter-annual climatic variability. High resolution satellite time series such as the one used in this study, are not yet available over multi-year periods but the advent of the Sentinel-2 mission will change this status quo and will provide more opportunities to address this issue.

Despite these limitations, we found that the model results are consistent with a semivariogram analysis of the DEM of the study area, since the semivariogram of the slope exhibits a break near 200 m. This result, in agreement with previous studies [START_REF] Essery | Scaling and parametrization of clear-sky solar radiation over complex topography[END_REF], suggests that a similar semivariogram analysis could be performed in other study areas as a guide to determine the model resolution. Further work is needed to confirm if this recommendation holds in regions with distinct climate and topography, but we expect that similar result could be derived in semi-arid mountain regions where the incoming solar radiation is a key factor in the snowpack dynamics.

Chapitre 4

Estimation du SWE dans le bassin versant de l'Ourika à partir de ré-analyses climatiques désagrégées

Ce chapitre se présente sous forme d'un article publié dans la revue Water : Baba, M. W., Gascoin, S., Jarlan, L., Simonneaux, V., and, Hanich., Variations of the Snow Water Equivalent in the Ourika Catchment (Morocco) over 2000[START_REF] Africanews | Morocco experiences heavy snowfall, social media reacts[END_REF]Using Downscaled MERRA-2 Data. Water 2018, 10(9), 1120, http: // doi. org/ 10. 3390/ w10091120

Motivations et résumé de l'étude

Depuis les années 2000, plusieurs études ont permis d'observer ou de modéliser le manteau neigeux dans le Haut Atlas marocain [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF][START_REF] Boudhar | Analysis of snow cover dynamics in the Moroccan High Atlas using SPOT-VEGETATION data[END_REF][START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. Or, jusqu'à présent, aucune étude n'a visé à estimer d'une façon distribuée l'équivalent en eau dans le Haut Atlas marocain sur une période pluriannuelle. De plus, la plupart des études se sont concentrées sur le bassin de la Rheraya car ce bassin est le site d'étude privilégié dans le cadre du programme SUDMED, de l'observatoire du Tensift et du LMI TREMA. Pourtant, la contribution hydrologique de l'oued Rheraya est inférieure à celle de l'oued Ourika (Fig 2 .4).

Le but de cette étude est de proposer la première estimation de l'état du manteau neigeux depuis 2000 jusqu'au printemps de 2018 dans le bassin de l'Ourika (Fig. 2.1). L'estimation est faite à l'aide du modèle du manteau neigeux SnowModel qui est alimenté par des données de réanalyses météorologiques MERRA-2. La résolution spatiale des simulations a été fixée à 250 m, d'après les résultats du chapitre 3 cela permet d'avoir une représentation correcte des flux radiatifs notamment avec un temps de calcul raisonnable. Cela nous a permis d'étendre la période de l'étude sur toute la période couverte par les données MODIS (18 ans), afin de caractériser la variabilité inter-annuelle du manteau neigeux sur une période assez longue tout en ayant un jeu de données indépendant pour évaluer la qualité des simulations. L'évaluation des simulations a été effectuée avec des données de hauteur de neige observées au sommet d'Oukaimeden (3230 m.a.s.l), les données de surfaces enneigées MODIS, et des données hydro-météorologiques in situ (précipitation et débit).

Les principales conclusions tirées de cette étude sont :

-Les simulations montrent une grande variabilité inter-annuelle du SWE, soit en terme de quantité de l'équivalent en eau, de surfaces enneigées et aussi de durées de l'enneigement. -La comparaison aux images de couverture neigeuse de MODIS pour la majorité des années a montré que le modèle reproduit généralement bien la distribution de la couverture neigeuse. Les précipitations sont moins bien reproduites probablement car les pluviomètres sont situés à plus basse altitude que les surfaces enneigées. De plus il semble que la qualité des simulations est meilleure en fin de période. -La qualité des simulations est correcte compte-tenu de la résolution spatiale des données MERRA-2 en entrée. Ceci suggère que ces données peuvent combler le vide causé par le manque de stations météorologiques dans le Haut Atlas si elles sont désagrégées. Des résultats similaires ont été obtenus dans d'autres massifs [START_REF] Cortés | Snow process estimation over the extratropical andes using a data assimilation framework integrating MERRA data and Landsat imagery[END_REF][START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF]. -Les données de réanalyses météorologiques MERRA-2 permettent de suivre l'enneigement avec un délai de 2 semaines.

Pour conclure, le modèle distribué de neige (SnowModel) alimenté par les données MERRA-2 est une solution intéressante pour fournir une première estimation de l'état du manteau neigeux. Toutefois, des biais importants sont surélevés et sont probablement dus aux erreurs de forçages météorologiques (principalement, les précipitations et la température). Ceci nous conduit dans le chapitre suivant à proposer un schéma d'assimilation de données qui vise à réduire ce biais.

Article

Abstract

The Ourika River is an important tributary of the Tensift River in the water-stressed region of Marrakesh (Morocco). The Ourika river flow is dominated by the snow melt contribution from the High Atlas mountains. Despite its importance in terms of water resources, the snow water equivalent (SWE) is poorly monitored in the Ourika catchment. Here, we used MERRA-2 data to run a distributed energy-balance snowpack model (Snow-Model) over 2000-2018. MERRA-2 data were downscaled to 250 m spatial resolution using a digital elevation model. The model outputs were compared to in situ measurements of snow depth, precipitation, river flow and remote sensing observations of the snow cover area from MODIS. The results indicate that the model provides an overall acceptable representation of the snow cover dynamics given the coarse resolution of the MERRA-2 forcing. Then, we used the model output to analyze the spatio-temporal variations of the SWE in the Ourika catchment for the first time. We suggest that MERRA-2 data, which are routinely available with a delay of a few weeks, can provide valuable information to monitor the snow resource in high mountain areas without in situ measurements.

Introduction

The snowpack is an important resource in many semi-arid and Mediterranean catchments, mostly because it provides meltwater runoff for irrigation during the crop growing season. Snow melt is also important to replenish reservoirs before the dry summer season and contribute to the natural recharge of the groundwater in the lowlands [START_REF] Viviroli | Mountains of the world, water towers for humanity : Typology, mapping, and global significance[END_REF][START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF].

The Tensift basin is a major basin in Morocco covering 20450 km 2 around the city of Marrakesh from the High Atlas mountains to the Atlantic coast. In this semi-arid region of Morocco, the precipitations are scarce while the evaporative demand is high [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. In the center of the basin, the Haouz plain (6000 km 2 ) is extensively cultivated and many farmers rely on irrigation to grow their crops. In Marrakesh, where the urban population is approximately one million inhabitants and rapidly growing (+2.5 % per year) (RGPH, 2014), the mean annual precipitation is only 250 mm and is concentrated between November and April [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. By contrast, due to the orographic enhancement of the precipitation, the mean annual precipitation measured at Oukaimeden-CAF (2650 m.a.s.l) in the High Atlas region of the Tensift basin, ranged between 300 mm and 900 mm over 1989[START_REF] Daget | Revue des méthodes d'assimilation[END_REF][START_REF] Simonneaux | Modelling runoff in the Rheraya Catchment (High Atlas, Morocco) using the simple daily model GR4J. trends over the last decades[END_REF]. This precipitation pattern drives the hydrology of the Tensift basin : most of the streamflow in the Tensift river originates from five main tributaries which have their headwaters in the northern slope of the High Atlas mountains [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]. These tributaries are (from east to west) the Nfis, Rheraya, Ourika, Zat and R'Dat rivers (Table 4. 1,Fig. 4.1). The estimated contribution of snow melt to streamflow in these catchments over the period 2002 to 2005 ranged between 2 % for the Rda't catchment in 2004 to 51 % in 2003 for the Ourika catchment [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]. Previous studies about the hydrological processes in the High Atlas have focused on the Rheraya catchment because this catchment was selected as a pilot site in the SudMed programme and LMI TREMA [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. These studies highlighted the importance of snow melt to the hydrological regime of the Rheraya river [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Richard | Evaluation d'un modèle pluie-debit sur le bassin versant de la Rheraya dans le Haut Atlas marocain[END_REF] and the large interannual variability of the snow accumulation over the last decades [START_REF] Simonneaux | Modelling runoff in the Rheraya Catchment (High Atlas, Morocco) using the simple daily model GR4J. trends over the last decades[END_REF]. In contrast, the snow regime in the Ourika catchment is largely undocumented, although this catchment arguably has a greater importance in terms of water resources due to Catchment Area (km its larger size. In addition, discharge measurements indicate that the Ourika catchment has the highest mean specific discharge among these five catchments (Table 4.1), which highlights its significance as a runoff production area.

2 ) Z 0 (m) z (m) Q (Mm3) Q S (
Previous studies have analyzed the interannual variability of the snow cover area (SCA) in the Ourika catchment (among other catchments) using MODIS data [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF], but the SCA is not a direct indicator of the water resources that is brought by the snowpack. The main variables of interest for water managers are the snow water equivalent (SWE) and snow melt.

There are different approaches to simulate SWE and snow melt in mountain catchments. Degree-day models are practical to implement since they require only the daily precipitation and temperature data [START_REF] Hock | Temperature index melt modelling in mountain areas[END_REF]. On other hand, physically-based distributed snowpack models solve the snowpack energy balance equation [START_REF] Marks | A spatially distributed energy balance snowmelt model for application in mountain basins[END_REF]. They offer a more realistic representation of the snow cover evolution especially in Mediterranean semi-arid areas where the solar radiation is a key forcing of the snowpack energy budget and the sublimation can represent a significant part of the ablation [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF][START_REF] López-Moreno | Different sensitivities of snowpacks to warming in Mediterranean climate mountain areas[END_REF][START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF]. Unlike degree-day models, energy balance models do not require the calibration of the model parameters, but they require more meteorological forcing data (precipitation, air temperature, air humidity, solar shortwave radiation, atmospheric longwave radiation, wind speed, air pressure). In many mountain regions like the High Atlas, these measurements are generally not available, and there is no direct way to measure the SWE in mountains by satellite remote sensing to this date [START_REF] Dozier | Estimating the spatial distribution of snow water equivalent in the world's mountains[END_REF], therefore a solution is to use atmospheric model outputs as forcing to run a snowpack model [START_REF] Brun | Simulation of northern eurasian local snow depth, mass, and density using a detailed snowpack model and meteorological reanalyses[END_REF][START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF][START_REF] Alonso-González | Daily gridded datasets of snow depth and snow water equivalent for the Iberian Peninsula from 1980 to 2014[END_REF]. In particular, climate model reanalysis like NASA modern-era retrospective analysis for research and applications data (MERRA-2) provide an optimal reconstruction of the meteorological variables at the global scale by assimilation of available in-situ and satellite observations into a climate model. Previous studies showed that MERRA-2 can provide a realistic representation of the precipitation for snow and hydrological studies [START_REF] Gelaro | The modern-era retrospective analysis for research and applications, version 2 (MERRA-2)[END_REF]. In Africa the MERRA-2 precipitation correction algorithm uses a specific configuration, which has been introduced to improve deficiencies in MERRA-Land (Reichle et al., 2017a). This configuration has improved the simulation of the terrestrial water storage in Africa at coarse resolution as observed by GRACE, but MERRA-2 skills were not specifically evaluated in North Africa [START_REF] Gelaro | The modern-era retrospective analysis for research and applications, version 2 (MERRA-2)[END_REF].

The objective of this study is to evaluate the combination of MERRA-2 data with a distributed snowpack model to simulate and analyze the spatio-temporal variability of the snow water equivalent (SWE) in the Ourika catchment from 01 November 2000 to 01 April 2018. The study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) was chosen because it corresponds to the observation period of the MODIS sensor, which we use as validation data to evaluate our SWE simulation.

We use SnowModel (Liston and Elder, 2006a), a distributed energy-balance model to simulate the snowpack evolution. A key asset of SnowModel is its MicroMet subroutine, which spatially interpolates the meteorological forcing over the domain using the digital elevation model. In addition, SnowModel accounts for the effect of the terrain slope and aspect on the incoming solar radiation, which can be significant in this semi-arid mountainous region. The MERRA-2 forcing were downscaled to 250 m using a digital elevation model to account for the effect of the topography on the temperature, humidity, wind, precipitation and solar radiation [START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF].

Finally, we take advantage of the availability of the MERRA-2 forcing in near real-time until April 2018 to analyze the snow condition during the winter 2017-2018. According to local newspaper, this winter was characterized by heavy snowfall events in the High Atlas above 1000 m in the middle of January and the beginning of February [START_REF] Africanews | Morocco experiences heavy snowfall, social media reacts[END_REF]. A rare snowfall event in the lower elevation areas also made the headlines in various media outlets (alj, 2018).

Data and method

Study area

The simulation domain covers the Ourika catchment at Aghbalou (982 m). The catchment is centered at 31 o 10N, 7 o 40W between the semi-arid Haouz plain and the northern slope of the High Atlas (Fig. 4.1) .

The Ourika river provides water resources to the rural communities living in the valleys of the High-Atlas (e.g. Setti-Fatma, Ourika). Meltwater is used to irrigate croplands and gardens (e.g. Paradise of saffron, Ourika Organic Aromatic Garden) and feeds many irrigation systems in the High-Atlas foothills and the Haouz plain until its confluence with the Tensift river near Marrakech.

The catchment elevation ranges between 974 m and 4001 m.a.s.l (Jebel n'Tarourt). The fraction of the catchment area above 3000 m.a.s.l is 25 %. Winter precipitation mainly falls as snow in these high-elevation areas, but snow is also present above 2000 m more intermittently [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF].

The catchment is sparsely vegetated, except for the valleys bottoms. The catchment geology favors the rapid runoff of snow melt and rainfall since two thirds of the catchment is occupied by Precambrian magmatic rocks characterized by low infiltration rate [START_REF] Juillerat | Etude hydrologique comparative de 5 bassins versants du Haut Atlas marocain[END_REF].

The Ourika river regime at Aghbalou is dominated by the snow melt contribution with a maximum discharge in April (Fig. 4.2). In average, 60 % of the annual discharge occurs between March and May (Agence Japonaise de coopération internationale, 2008) (Fig. 4.2). However, the hydrological regime is also influenced by the high annual variability in the snowfall amount and melt-out date [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF][START_REF] Zkhiri | Regional flood frequency analysis in the High Atlas mountainous catchments of Morocco[END_REF]. to the daily time step for this study. There are many gaps in the measurements but they still enable to evaluate the model. The available data covers six snow seasons (2004-2005, 2005-2006, 2008-2009, 2010-2011, 2013-2014 and 2014-2015).

The meteorological forcing to run SnowModel were derived from the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2), available at 0.5 o × 0.625 o of spatial resolution. Each new month of MERRA-2 reanalysis is available approximately between the 15th and the 20th of the next month. The extracted variable are listed in Tab. 4.2. Each MERRA-2 pixel is characterized by its latitude and longitude and surface geopotential height (named as PHIS in the data file), which were converted to (x, y, z) coordinates in the WGS-84 UTM-29N projection system (see subsection 4.2.3).

A digital elevation model (DEM) is required to downscale the meteorological forcings 3). Here we used the Shuttle Radar Topography Mission Void Filled (SRTM-Void Filled Version 4.1) [START_REF] Jarvis | Hole-filled SRTM for the globe version 4. available from the CGIAR-CSI SRTM 90m Database[END_REF]. It was resampled with the bicubic method to 250 m spatial resolution to limit to computational cost of the simulation. This resolution was found to be a good trade-off between the computational cost and the accuracy of the simulation, because it allows capturing the main snow cover heterogeneities caused by the different radiative energy at the hillslope scale [START_REF] Baba | Contribution of high resolution remote sensing data to the modeling of the snow cover the in Atlas Mountains[END_REF].

The catchment land cover was obtained from ESA-CCI classification map (Fig. 4.3) [START_REF] Bontemps | Consistent global land cover maps for climate modelling communities : current achievements of the ESA's land cover CCI[END_REF]. The majority of the region (especially the high elevated areas) is classified as bare soil.

We used 18 years of daily MODIS data at 500 m of spatial resolution to validate the simulation. These data provide the snow cover fraction (SCF) within each MODIS pixel. The cloud pixels were interpolated using the algorithm of [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF], which was specifically developed for the High Atlas. The accuracy of the algorithm was evaluated by comparison to snow depth measurements and high-resolution snow mask derived from the Formosat-2 observations over the Rheraya catchment. The results showed that the algorithm allowed a reduction of the cloud coverage to less than 1% with an overall accuracy of 89% [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. The MODIS data were produced from the collection 5 MOD10A1 snow product. However, the collection 5 has been discontinued since 01 January 2017 therefore the entire time series was reprocessed using collection 6 MOD10A1. Since the collection 6 does not provide anymore the SCF but only the normalized difference snow index (NDSI) of snow-covered pixels, we applied the same NDSI-SCF relationship that was used in the collection 5 [START_REF] Hall | Algorithm theoretical basis document (ATBD) for the modis snow and sea ice-mapping algorithms[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF] : SCF = -0.01 + 1.45 × NDSI (4.1)

Then, we applied the cloud-removal algorithm to generate a continuous time series of the SCF from 01 September 2000 to 01 April 2018 in our study area that is fully consistent with the data presented by [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF].

We used monthly river discharge data at the Aghbalou river gauge (982 m, Fig. 4.2) from the Tensift basin agency from 2000 to 2017 during the melt season (January to May). The daily discharge data were only available for 7 years in this period. We used the daily data to compute the average discharge over a 10-day time step during the melt season.

We used daily precipitation measurements from three automatic weather stations (AWS) in the Ourika catchment (Fig.

Methodology

To simulate the SWE, the melt and the discharge in the Ourika catchment, we fed SnowModel, a distributed energy balance snowpack evolution model (Liston and Elder, 2006a) with MERRA-2 reanalysis. The SnowModel has previously been tested with success in the Arctic, Antarctic, and high mountain regions, comparing simulated snow accumulation, distribution, and ablation processes with observations [START_REF] Hiemstra | Snow redistribution by wind and interactions with vegetation at upper treeline in the Medicine Bow Mountains, Wyoming, USA[END_REF][START_REF] Boulaajoul | Modelling the snow distribution at two high arctic sites at Svalbard, Norway, and at an alpine site in central Norway[END_REF][START_REF] Hiemstra | Observing, modelling, and validating snow redistribution by wind in a wyoming upper treeline landscape[END_REF][START_REF] Mernild | Climatic control on river discharge simulations, Zackenberg River drainage basin, northeast Greenland[END_REF][START_REF] Mernild | Runoff and massbalance simulations from the greenland ice sheet at Kangerlussuaq (Søndre Strømfjord) in a 30-year perspective, 1979-2008[END_REF][START_REF] Suzuki | Impact of land-use changes on snow in a forested region with heavy snowfall in Hokkaido, Japan[END_REF][START_REF] Mernild | Greenland freshwater runoff. part ii : Distribution and trends, 1960-2010[END_REF][START_REF] Gascoin | Wind effects on snow cover in Pascua-Lama, Dry Andes of Chile[END_REF][START_REF] Szczypta | Impact of climate and land cover changes on snow cover in a small Pyrenean catchment[END_REF][START_REF] Suzuki | Variations of winter surface net shortwave radiation caused by land-use change in northern Hokkaido, Japan[END_REF]. SnowModel is composed of four sub-models : MicroMet (meteorological data interpolation), EnBal (energy balance computations) (Liston and Elder, 2006a), SnowPack (snow depth evolution) and SnowTran-3D (wind transport) [START_REF] Liston | Simulating complex snow distributions in windy environments using SnowTran-3D[END_REF]. Given the lack of in situ meteorological observations, SnowTran-3D was not activated. We set all the SnowModel parameters to default values because we have a limited observations dataset in the study area. By using only default parameters we can evaluate the accuracy of the model if it should be transposed to other catchments in the High Atlas without calibration.

MicroMet is the key interface between the MERRA-2 data and the snowpack model, since it enables to distribute all meteorological data that are required to run EnBal and SnowPack sub-models over the study area at the resolution of the DEM. These variables are precipitation, air temperature, relative humidity, wind speed, incoming solar radiation and incoming longwave radiation at the hourly time step. MicroMet was originally designed to interpolate station data onto a regular grid, but it can also be applied to downscale large-scale, surface-level forcing from climate model or reanalysis such a MERRA-2 [START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF]. In this case, each MERRA-2 grid cell intersecting the study domain is considered as a virtual station, which is located at the center of the cell.

To distribute the meteorological data from (virtual) station data, MicroMet performs a horizontal interpolation using the Barnes objective analysis scheme, which applies a Gaussian distance-dependent weighting function to interpolate the station data to the model grid resolution. Topography-dependent functions are applied to account for the effect of elevation, slope and aspect on the gridded meteorological forcing. In the case of air temperature, a linear lapse rate correction is used :

T x = T stn -τ (Z x -Z stn ), (4.2)
where T stn ( o C) is the air temperature at elevation Z stn (m), T x ( o C) refers to the temperature at a given elevation Z x (m), and τ is default monthly lapse rate ( o C.m -1 ) [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. A non-linear function is used for the precipitation :

P x = P 0 × [1 + χ(Z x -Z 0 )] [1 -χ(Z x -Z 0 )] , ( 4.3)
where P x (mm) is the precipitation at Z x , P 0 (mm) is the interpolated precipitation, Z 0 (m) is the interpolated station elevation and χ (km -1 ) is the monthly precipitation correction factor which we set to default values [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. The distributed relative humidity is derived from similar rules, while the distributed wind speed is generated by accounting for the deflection of the wind field by the terrain. The longwave and shortwave radiations are interpolated using the Barnes scheme, but the shortwave radiation is modified according to the slope and orientation. For further details the reader can refer to (Liston and Elder, 2006a).

SnowModel takes as input the coordinates of the stations. However, a pre-processing of MERRA-2 variables was developed to make them compatible to SnowModel input. LON and LAT were projected to x, y coordinates in the WGS-84 UTM 29N projection system. The elevation of the MERRA-2 grid cell was computed from the geopotential height (m 2 /s 2 ) by dividing it by the Earth's gravity (m/s 2 ).

The relative humidity RH was retrieved from QV2M as follows : RH = Q2M 0.62e s (T ) (4.4) Where e s (T ) is the water vapor saturation pressure, which is a function of temperature and approximated by Magnus-Tetens formula e s (T ) = 0.61exp( 17.625T T +243.04 ). Eastward (u) and northward (v) wind components were converted to wind speed (WS) and wind direction (WD) using :

WS = √ u 2 + v 2 (4.5) WD = arctan(v/u) (4.6)
The validation of the model outputs is challenging due to the lack of SWE measurements in Morocco and the difficulty of getting it otherwise (e.g the cost of LIDAR campaigns [START_REF] Kirchner | LiDAR measurement of seasonal snow accumulation along an elevation gradient in the southern Sierra Nevada, California[END_REF][START_REF] Nadal-Romero | The application of terrestrial laser scanner and SfM photogrammetry in measuring erosion and deposition processes in two opposite slopes in a humid badlands area (central Spanish Pyrenees)[END_REF] is very high, and the use of unmanned aerial vehicles (UAV) to estimate the snow depth [START_REF] Boesch | Comparison of digital surface models for snow depth mapping with UAV and aerial cameras[END_REF] is limited by its spatial coverage on one hand and the complicated legal procedure to use it in Morocco on other hand. To deal with this, we have combined different hydrometeorological datasets and remote sensing observations. The simulated snow depth (SD) was extracted at the Oukaimeden station (3200 m.a.s.l) to assess the model performance using the Root Mean Square Error (RMSE). We also evaluate the cumulative precipitation over the snow season (from 1 st December to 1 st April from 2000 to 2009 by using in-situ precipitation. The model was further evaluated at the catchment scale using the daily snow cover area of the Ourika catchment (SCA) from MODIS data (Fig. 4.4). The simulated SCA was computed from the simulated SWE by computing, for every day of the simulation period, the number of grid cells where the simulated SWE exceeds a fixed threshold of 20 mm [START_REF] Gascoin | A snow cover climatology for the Pyrenees from MODIS snow products[END_REF]. Then, the snow covered area of the catchment was computed and compared to the MODIS data using the coefficient of determination (R 2 ) for each hydrological year from 01 September to 31 May. We excluded the period from 01 June to 31 August, because the snow is absent during these months, which could artificially increase the R 2 . We verified that varying a SWE-SCA threshold between 20 mm to 50 mm did not significantly change the results. We also used the Ourika discharge to evaluate the simulated runoff. SnowModel does not have a soil reservoir and only computes for each grid cell a "runoff", which is the total liquid water that leaves the snowpack or the rainfall if there is no snow cover. The liquid water in the snowpack comes from the snow melt or the rain-on-snow. It is released when the liquid water content exceeds the snowpack maximum storage capacity. Given the lack of hydrological routing in the model, we compared the modeled runoff with the observed discharge using monthly averages only.

Results

Validation

In general, the model simulated reasonably well the snow depth at Oukaimeden station (Fig. 4.5). The model failed to properly simulate the evolution of the snow depth in 2004-2005 because of inconsistent precipitation events. The accumulation was underestimated in 2013-2014, but otherwise there is no clear tendency to overestimate or underestimate the snow depth at this point (Tab. 4.3). The melting rates were correctly reproduced in 2005-2006 and 2008-2009. The comparison with MODIS data shows that the model was able to produce a reasonable simulation of the snow cover area in the catchment (Fig. 4.6). We found a R 2 greater than 0.70 for the majority of the hydrological years (Tab. 4.4).

However, the SCA was poorly simulated in 2000-2001, 2002-2003, 2005-2006 and 2006-2007. We note that the model performance tends to increase over the study period. In particular, all R 2 are higher than 0.77 between 2009-2010 and 2016-2017. The R 2 is slightly lower in 2017-2018 (R 2 =0.73) but the snow season was not completely over at the time of writing this study.

The simulated precipitation at Oukaimeden-CAF is poorly correlated with the observations (R 2 = 0.22), expect for the most recent months (Fig. 4.7). At Amenzal and Agouns stations, the simulated precipitation does not agree with the observations, and contrary to Oukaimeden-CAF, the model largely overestimates the precipitation in the recent months (Fig. 4.7). The model strongly overestimates the precipitation in 2002-2003 Table 4.4 -R 2 between the simulated and observed snow cover area (SCA) and discharge (Q) in the Ourika catchment for every hydrological year of the study period. The R 2 was computed at the daily timestep from 01 September to 31 May for SCA and at the monthly timestep from January to May for Q (NA : no measurements available for this year).

Years .9 show the comparison between the simulated runoff and the observed discharge at the outlet of the Ourika catchment. The mean R 2 between the modeled runoff and the observed discharge at the monthly time step from 2000 to 2017 is 0.55 ; 10 on 15 years have a R 2 value higher than 0.5 (Fig. 4.8). However there remain large discrepancies between both datasets, in particular the low R 2 in 2000-2001, 2001-2002, 2003-2004 and 2004-2005 are consistent with the lower R 2 which were obtained for this year for the SCA. However, this is not always the case, as for instance a low R 2 in discharge was obtained in 2012-2013 despite a correct simulation of the SCA for this hydrological year (Fig. 4.6). Fig. 4.9 suggests that at the finer 10-day time step, the model can capture peak discharges, however, the discharge regime is generally poorly represented except in 2012 and 2015.

Variations of the simulated SWE

The simulated SWE was averaged over the catchment area and plotted in Fig. 4.10. The catchment SWE typically ranged between 80 and 120 mm we during the winter (Fig. 4.10). However, this range was largely exceeded, (e.g. in 2006-2007 and 2014-2015) with values reaching 300 mm and 400 mm. Fig. 4.10 indicates that the melt season generally starts in March since the slopes of the SWE time series are generally negative after March. But again there is a large variability. Indeed,Fig. 4.11 (top panel) indicates that the annual maximum SWE,an indicator of the maximum water reserve in the catchment, varied between 290 mm we and 1100 mm we over the study period. The date of the peak SWE also exhibited a large variability, from early March to mid May (Fig. 4.11,middle panel). The date of the snowpack snow disappearance date has the same range of variability from mid March to late May (Fig. 4.11,lower panel). Using the Mann-Kendall test, we did not find any significant trend in these indicators at the 5 % confidence level.

In order to visualize the spatial distribution of the simulated SWE within the catchment we extracted the April 01 SWE for every year (Fig. 4.12). The SWE on April 01 is commonly used as an indicator of the annual SWE before the melt season in the Northern Hemisphere [START_REF] Dewalle | Principles of snow hydrology[END_REF]. In contrast with the large temporal variability of the SWE, there is consistent pattern in the SWE spatial distribution at the scale of the catchment. Again, we did not find any significant trend at the resolution of the grid scale in the April 01 SWE. The model suggests that the SWE can reach 1000 mm we on April 01 in the highest part of the catchment during the wettest years (2005-2006, 2008-2009, 2014-2015). However, the SWE can be almost null in other years (Fig. 4.12).

In the beginning of the 2017-2018 snow season, the simulated SWE in the Ourika catchment was amongst the lowest values since 2000 (Fig. 4.10). However, after the large snowfall events in January and February, the SWE exceeded the mean value of the previous 17 years. Despite this high SWE value in mid February, the snowpack returned to a record low value on April 01 after one month of intense melt. The SWE was very low even at high elevation on April 01 (Fig. 4.10).

Discussion

The evaluations of the simulation at the point-scale and at the catchment scale indicate that the model performs rather well given the coarse resolution of the MERRA-2 input data. The point-scale comparison with the observed snow depth at Oukaimeden gives correct results except for 2004-2005. A single station has a limited spatial representativeness at the scale of the catchment, but it is the only direct measurements of the snow depth in the area and it provides a robust reference to evaluate the temporal variations of the simulated snow depth. The model outputs do not agree as well with the precipitations. The simulated precipitations are more consistent with observations at the high elevation station of Oukaimeden than with observations at Agouns and Amenzal stations. We think this is due to the fact that Agouns and Amenzal stations are located in valleys bottoms, while the Oukaimeden-CAF station is located in an open area. This station is probably less influenced by the local scale effects. Local atmospheric effects like orographic precipitation shadows are not resolved by a large scale atmospheric model such as MERRA-2, and are not represented in MicroMet either, since it relies on a standard precipitation lapse rate to distribute the precipitation [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF].

The validation of the simulated gridded precipitation with point-scale precipitation measurements remains challenging because the precipitation gauges are influenced by local wind turbulences, which can cause snow under-catch [START_REF] Fassnacht | Estimating Alter-shielded gauge snowfall undercatch, snowpack sublimation, and blowing snow transport at six sites in the coterminous USA[END_REF]. MODIS data suggest that the model performances at the catchment scale are better than what could be inferred from the precipitation measurements alone (Fig. 4.6). The simulated runoff is also relatively consistent with the river discharge for the majority of years. This indicates that the mean SWE in the catchment is generally well simulated.

MODIS data enable to better understand the causes of the model errors during those snow seasons, which exhibit the largest errors in the runoff simulation.

-In 2000-2001, the model produces high runoff values while the discharge was low.

In the same period, large variations in the SCA were simulated, while MODIS data indicate that the SCA was close to zero (Fig. 4.6). -In 2003-In -2004, the simulated SCA is well correlated to the MODIS SCA during the melt season, although the model overestimates the runoff from May to April. This suggests that there is an excess of liquid precipitation in the model forcing. -In 2008-2009, the runoff simulation is completely anti-correlated with the discharge measurements. Given that the snow depth at Oukaimeden and the SCA over the catchment are properly simulated (Fig. 4.5,4.6), we conclude that the simulated rise in runoff in spring is due to erroneous rainfall events in MERRA-2.

We have used the runoff as an indicator of the discharge since the model does not simulate the water routing and storage in the river network, subsurface flow and evaporation. Additional work is needed to account for the hydrological routing of the snow melt and rain runoff within the catchment to explore the model performance more deeply. However, modest performances with a calibrated hydrological model were reported by (Y. et al., 2018) in the neighboring catchment of the Rheraya river, although the Rheraya catchment is better monitored with several weather stations. The validation of discharge in the High Atlas catchments is difficult, ought to the frequent occurrence of flash floods, which modify the river bed and cast doubt on the rating curve accuracy.

The simulation yielded better results over the last decade, which suggests that the MERRA-2 reanalysis was improved over the Ourika region, by incorporating additional measurements in the recent years (Gelaro et al., 2017a). This confirms that this SWE simulation is not adapted for trend analysis, since the number of observations in the MERRA-2 assimilation scheme is not guaranteed to be constant. Therefore we cannot firmly conclude that there is no trend in the snow cover in the Ourika catchment from our analyses. In addition, any trend might be extremely difficult to detect given the large interannual variability of the snow cover in the Ourika catchment.

The model outputs highlight the large interannual variability of the SWE in the Ourika catchment, in agreement with previous studies in the High Atlas [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF][START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF][START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. This variability is typical of semi-arid catchments, and is much higher than what is observed in other Mediterranean mountains like the Pyrenees [START_REF] Gascoin | A snow cover climatology for the Pyrenees from MODIS snow products[END_REF] or Mount Lebanon [START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF].

Regarding the snow conditions in 2017-2018, the model indicates that this snow season was characterized by a series of extreme events. The winter started like a snow drought, with a long period without snow accumulation in November and December. The SWE rose drastically in January and February after two major snowfall events but melted rapidly afterwards (Fig. 4.10). The analysis of the air temperature over the Ourika catchment suggests that November 2017 was the hottest November since 2000, while February 2018 was the third coldest February since 2001. This low temperature in February 2018 explains the exceptional presence of snow in the low elevation areas. The rapid melt thereafter can be explained by a succession of hot days ; for instance, in March 2018, the spatiallyaveraged maximum temperature was 8.5 o , while the average for the period 2000-2018 is 8.0 o with a standard deviation of 1.2 o .

Summary and conclusions

We showed the application of a distributed snowpack model in the snow-dominated Ourika catchment in the High Atlas. We focused on the Ourika catchment because there is a lack of knowledge regarding the snow dynamics in this catchment although snow melt is an essential water resource to sustain irrigation in the semi-arid Haouz plain. To our best knowledge, this study is the first to evaluate and analyze spatially-distributed simulation of the SWE in the High Atlas mountains for a period covering almost the last two decades.

To cope with the lack of weather station measurements in the High Atlas, we used the MERRA-2 reanalysis and the MicroMet model to generate a distributed meteorological dataset, taking into account the main effects of topography on the meteorological variables. The model outputs were compared to a snow depth record and the time series of the catchment snow cover area from MODIS. The results suggests that the model gave a reasonable estimation of the snowpack evolution over the period 2000-2018, despite the coarse resolution of the MERRA-2 forcing. We have used SnowModel in its default configuration without calibration, which suggests that similar performances could be expected if the model is applied to other ungauged catchments in the High Atlas.

The simulated SWE exhibited a significant interannual variability, both in terms of amplitude (maximum amount of water stored in the snowpack) and duration of the snow season. The model also gave clues about the snow conditions during the hydrological year 2017-2018, which was characterized by different extreme meteorological events.

The model outputs, however, are still subject to large uncertainties.These results should be taken with caution since we could not use direct SWE measurements to validate the model, but only indirect observations of the snow cover dynamics. We suggest to conduct SWE measurement campaigns in the following years to better validate the SWE simulation over the catchment. Further work is required to evaluate the simulated SWE against discharge measurements at finer temporal timescales. The biases in the MERRA-2 forcing could also be reduced by assimilating snow cover maps retrieved from remote sensing data [START_REF] Margulis | A Landsat-era Sierra Nevada snow reanalysis (1985-2015)[END_REF][START_REF] Girotto | Examining spatial and temporal variability in snow water equivalent using a 27 year reanalysis : Kern River watershed, Sierra Nevada[END_REF]. Satellite-derived snow cover maps are usually available with a shorter delay than MERRA-2 (less than a week). Therefore it should be possible to develop a near-real time data assimilation scheme to produce more accurate SWE maps. This is the focus of our ongoing work.

Yet, we think that the simple approach shown here could already be useful for water managers. The availability of the MERRA-2 forcing with a short delay (about 1 month) now enables to produce interesting diagnostics regarding the state of the snowpack, over the course of the snow season, without ground measurements. These data could help assess the snow water resource in the Ourika catchment, and in other catchments in the High Atlas, to optimize the planning of the irrigation during the spring.

Introduction

Motivations du travail

Nous avons vu dans le chapitre précédent (Chap. 4) les apports et les limitations de la modélisation du manteau neigeux en utilisant un modèle distribué de l'évolution du manteau neigeux (SnowModel), alimenté par des forçages météorologiques globaux MERRA-2. Malgré de bonnes performances, les forçages sont parfois fortement biaisés ce qui peut conduire à des erreurs importantes sur l'estimation du SWE en particulier. Le but de ce chapitre est de développer un schéma d'assimilation de surfaces enneigées obtenues par télédétection spatiale à haute résolution spatiale et temporelle dans le modèle distribué de neige pour réduire le biais des forçages météorologiques et mieux estimer le SWE dans le bassin de la Rheraya.

Dans un premier temps, nous présenterons un aperçu général de l'usage de l'assimilation de données à travers le temps, puis en se concentrant sur les applications en nivologie. Ensuite, nous exposerons les ingrédients nécessaires pour concevoir et implémenter un filtre particulaire. Ce chapitre se termine sous forme d'un article soumis à Remote Sensing qui décrit le schéma d'assimilation proposé pour le suivi du manteau neigeux dans le Haut Atlas et montre son évaluation.

Une courte histoire de l'assimilation de données

Définition L'assimilation de données peut être définie comme un procédé qui consiste à combiner les observations de l'état d'un système donné avec les estimations de l'état obtenues par un modèle, de façon à avoir une estimation postérieure qui est plus proche de la réalité.

Imaginez vous dans un match de football ! Lors d'un coup de pied arrêté, le gardien de but prédit où le joueur va tirer le ballon (modèle). Au fur et à mesure, le gardien voit le ballon prendre sa trajectoire (observations) si bien qu'il modifie sa prédiction pour qu'elle s'approche de l'état réel (analyse). C'est cela l'assimilation de données : nous possédons un modèle et des observations ! Et l'assimilation de données permet de combiner ces deux ressources pour obtenir une meilleure estimation de la réalité.

C'est quoi l'assimilation de données selon l'auteur du manuscrit ? Histoire L'assimilation de données est souvent associée au développement des ordinateurs et des calculateurs scientifiques. En réalité, cette technique remonte au XVIII e siècle, où elle fut utilisée par Tobias Mayer (1723-1762), un astronome allemand [START_REF] Daget | Revue des méthodes d'assimilation[END_REF]. Il se servit de cette méthode pour réduire les erreurs des prédictions des mouvements de la Lune. Toujours au XVIII e siècle, le célèbre mathématicien et astronome allemand Karl Friedrich Gauss (1777[START_REF] Kyle | The status of the Nimbus-7 earth-radiation-budget data set[END_REF] parvint à prédire avec succès le mouvement de l'orbite de Cérès, sans passer par la résolution des équations non linéaires de Kepler décrivant les mouvements planétaires. Il a utilisé un ensemble limité d'observations incertaines pour calculer l'intervalle de l'orbite de l'astéroïde. Gauss introduisit donc l'idée que le modèle peut servir à l'estimation d'un état antérieur (prior en anglais), tandis que les observations corrigent cette estimation pour la rapprocher de la réalité, dite postérieure (posterior en anglais). Ce travail permit de jeter les bases de la méthode des moindres carrés, qui a été introduite par Adrien Marie Legendre (1752-1834). Par la suite et jusqu'au début du XX e siècle, les méthodes d'assimilations servirent principalement dans le domaine de l'astronomie [START_REF] Daget | Revue des méthodes d'assimilation[END_REF]. Au début du XX e siècle, Ronald Aylmer Fisher (1890-1962), un biologiste et statisticien britannique introduisit de nombreux concepts qui donnèrent naissance aux statistiques modernes. L'un de ces concepts est le maximum de vraisemblance (maximum likelihood, en anglais) qui est l'un des socles de la théorie de l'estimation. Après la deuxième guerre mondiale, l'essor de l'informatique permit le développement des méthodes statistiques. Ainsi, la théorie de l'estimation atteignit sa maturité dans les années 1960 [START_REF] Daget | Revue des méthodes d'assimilation[END_REF]. Durant cette même période, Rudolf Kalman (1930Kalman ( ,2016) ) introduisit les filtres séquentiels dont l'efficacité fut démontrée lorsque Kalman parvint à estimer l'évolution du module d'Apollo par rapport à la Lune. La figure 5.1 illustre comment l'usage du filtre de Kalman et des observations ont permis d'obtenir une estimation plus réaliste de la position du module d'Apollo (ce qu'on appelle l'analyse). Aujourd'hui, ce filtre est utilisé pour différentes technologies (géolocalisation, traitement d'images, téléguidage, etc.) et sciences. En particulier l'assimilation de données est couramment utilisée dans le domaine de la météorologie et des sciences de l'environnement. La première application remonte aux années 1950 lorsque George Cressman développa une approche d'interpolation appelée analyse objective ou analyse de Cressman permettant de rapprocher le modèle des observations (Fig. 5.2) [START_REF] Blayo | Introduction to data assimilation[END_REF]. En pratique, cette méthode possède plusieurs limitations 1) l'absence d'une méthode directe pour spécifier le poids optimal, 2) les erreurs des observations ne sont pas prises en compte, donc une assimilation de mesures erronées peut dégrader les prévisions. Malgré ces limitations, cette analyse est encore utilisée par sa simplicité et sa rapidité.

La disponibilité de nombreuses données météorologiques observées en temps réel a motivé le développement de nouvelles techniques d'assimilation de données plus performantes. Plus précisément [START_REF] Sasaki | Some basic formalisms in numerical variational analysis[END_REF] proposa les méthodes variationelles afin de prévoir la trajectoire des ouragans [START_REF] Daget | Revue des méthodes d'assimilation[END_REF]. Ces approches sont largement utilisées dans les centres météorologiques comme le centre européen pour les prévisions météorologiques à moyen terme ECMWF (en anglais : European Centre for Medium-Range Weather Forecasts) [START_REF] Derber | A reformulation of the background error covariance in the ECMWF global data assimilation system[END_REF][START_REF] Bonan | Assimilation de données pour l'initialisation et l'estimation de paramètres d'un modèle d'évolution de calotte polaire[END_REF]. L'approche variationelle peut s'appliquer en 3 dimensions (X,Y,Z) ou en 4 dimensions en ajoutant la dimension temporelle, elles sont respectivement nommées 3D Var et 4D Var [START_REF] Blayo | Advanced Data Assimilation for Geosciences : Lecture Notes of the Les Houches School of Physics : Special Issue[END_REF]. Ces techniques se basent sur la minimisation de la distance entre la trajectoire des observations et celle du modèle. La façon avec laquelle ces méthodes sont définies les rend souples et permet d'introduire facilement de nouveaux types observations dans le système d'assimilation

Utilisation de l'assimilation des données en hydrologie nivale

En hydrologie nivale, plusieurs méthodes d'assimilation de données sont utilisées. Chacune d'entre elles est caractérisée par des avantages et des limitations en fonction de son degré de complexité et du modèle de neige utilisé. L'une des approches les plus simples est l'insertion directe. [START_REF] Liston | A simple data assimilation system for complex snow distributions (SnowAssim)[END_REF] l'ont appliquée dans SnowModel en intégrant les observations du SWE afin d'ajuster les précipitations et la fonte. A travers les résultats obtenus, ils ont conclu que cette approche est efficace que pour des études très localisées. Néanmoins, si le nombre d'observations du SWE est limité, cette approche ne permettra pas de bien reproduire la distribution réelle de la neige. A ce jour, cette méthode est encore utilisée dans des centres de prévisions météorologiques qui se basent sur cette méthode d'assimilation pour produire des estimations des variables nivologiques. Ces dernières, sont par la suite utilisées pour lancer les nouvelles prévisions [START_REF] Charrois | Assimilation de réflectances satellitaires du domaine visible et proche infrarouge dans un modèle détaillé de manteau neigeux[END_REF]. D'autres approches plus sophistiquées, qui sont d'ailleurs coûteuses en temps d'exécution, sont développées et testées dans les zones montagneuses. Les filtres d'ensemble de Kal-man et le filtre particulaire sont parmi les approches les plus répandues [START_REF] Charrois | Assimilation de réflectances satellitaires du domaine visible et proche infrarouge dans un modèle détaillé de manteau neigeux[END_REF]. Durand et al. (2008a) ont développé une méthode basée sur l'EnKF (Ensemble Kalman Filter en anglais) alimentée par des données de surfaces enneigées afin de perturber les précipitations et le rayonnement absorbé pour reconstruire le SWE. Ils ont démontré que cette approche pourrait être très utile dans les régions où le manteau neigeux est intermittent. D'autres études visaient aussi à coupler les méthodes d'ensemble de Kalman avec les données de surfaces enneigées dans le but de mieux estimer le SWE [START_REF] Clark | Assimilation of snow covered area information into hydrologic and land-surface models[END_REF][START_REF] Andreadis | Assimilating remotely sensed snow observations into a macroscale hydrology model[END_REF][START_REF] Su | Enhancing the estimation of continental-scale snow water equivalent by assimilating MODIS snow cover with the ensemble kalman filter[END_REF][START_REF] Leisenring | Snow water equivalent prediction using Bayesian data assimilation methods[END_REF][START_REF] Zhang | Assimilation of MODIS snow cover through the data assimilation research testbed and the community land model version 4[END_REF]. Une conclusion unanime est retirée de ces études : Une bonne performance de ces approches vis-à-vis de la reproduction de l'état du manteau neigeux est observée. En dépit de cette performance, ces méthodes sont confrontées principalement aux problèmes de linéarité [START_REF] Moradkhani | Uncertainty assessment of hydrologic model states and parameters : Sequential data assimilation using the particle filter[END_REF][START_REF] Piazzi | A particle filter scheme for multivariate data assimilation into a point-scale snowpack model in an Alpine environment[END_REF]. Le filtre particulaire permet de contourner cette limitation. En effet, il est très efficace pour les modèles non linéaires qui distribuent les sorties d'une manière non gaussiennes [START_REF] Moradkhani | Uncertainty assessment of hydrologic model states and parameters : Sequential data assimilation using the particle filter[END_REF][START_REF] Van Leeuwen | Nonlinear data assimilation in geosciences : an extremely efficient particle filter[END_REF]. [START_REF] Leisenring | Snow water equivalent prediction using Bayesian data assimilation methods[END_REF] ont comparé les performances du filtre particulaire avec celle d'un filtre d'ensemble de Kalman (EnKF). En utilisant ces deux approches, ils ont assimilé des observations synthétiques du SWE pour corriger les erreurs de forçages météorologiques (température et précipitations) afin de mieux reproduire l'évolution du manteau neigeux. Les résultats de cette comparaison suggèrent que chacune de ces méthodes permet de bien simuler cette évolution. En outre, en utilisant la même taille de chaque ensemble de simulations, ils ont remarqué une supériorité du filtre particulaire lors de l'estimation du SWE. Dans le même contexte, [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF] ont développé une approche de type (Particle Batch Smoother, PBS) qui est alimentée par les fractions de surfaces enneigées dans chaque pixel de 90 m pour estimer le SWE sur une période de 27 ans dans la Sierra Nevada (États-Unis). Les simulations du SWE en appliquant ce filtre ont été confrontées aux simulations du SWE par la méthode EnBS2 . Les résultats de cette comparaison ont permis de conclure que l'erreur quadratique moyenne du SWE est inférieure lors de l'usage du filtre particulaire. En plus de ces comparaisons, plusieurs études se sont basées sur l'alimentation d'un filtre particulaire avec des données de surfaces enneigées obtenues par télédétection spatiale pour reproduire le SWE dans les zones montagneuses [START_REF] Thirel | Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter[END_REF][START_REF] Cortés | Snow process estimation over the extratropical andes using a data assimilation framework integrating MERRA data and Landsat imagery[END_REF]. En particulier, [START_REF] Cortés | Snow process estimation over the extratropical andes using a data assimilation framework integrating MERRA data and Landsat imagery[END_REF] ont assimilé les fractions de surfaces enneigées dérivées de Landsat et reproduit l'évolution du SWE dans différents bassins versants des Andes en Amérique du Sud entre 1985 et 2015. Ils ont utilisé comme forçages météorologiques les données de réanalyses MERRA 3 . Les résultats de leur étude soulignent une amélioration marquante du SWE après la phase de l'assimilation. A la lumière de ces résultats, [START_REF] Cortés | Snow process estimation over the extratropical andes using a data assimilation framework integrating MERRA data and Landsat imagery[END_REF] suggèrent que la combinaison des données de réanalyses météorologiques et des données de télédétection peut apporter un bénéfice significatif dans les régions montagneuses où les stations météorologiques sont rares.

The results shown in this study highlight the usefulness of incorporating globally available atmospheric reanalysis and remotely sensed data sets into hydrological modeling. In the context of snow data products, the information generated from integrated data sets could have a significant impact in data scarce regions where coarse resolution climatological products are the only long-term, distributed forcing data set available for forward modeling exercises.

Cortés et al. (2016)

Pour synthétiser, le filtre particulaire alimenté par les données de surfaces enneigées et un modèle distribué de neige, permet de reproduire l'évolution du SWE. De plus, ce filtre est indépendant de la linéarité du système [START_REF] Van Leeuwen | Nonlinear data assimilation in geosciences : an extremely efficient particle filter[END_REF]. Cette étude bibliographique nous a conduit à choisir le filtre particulaire afin d'estimer l'état du manteau neigeux dans le Haut Atlas marocain.

Avant de décrire son implémentation, nous présentons ci-dessous quelques notions de bases utiles pour la suite.

Théorie de l'assimilation de données

Notations Vecteur d'état

Noté x(t), ce vecteur contient des informations concernant l'état physique du système étudié à l'instant t. Il s'agit d'une discrétisation des variables continues du système (par ex. le SWE). L'objectif de l'assimilation de données est de mieux estimer ce vecteur d'état comme nous le verrons par la suite. Notons que lorsque nous l'étudierons avec un temps discrétisé (t k ), ce vecteur est noté x k .

Modèle dynamique (SnowModel)

Pour simuler l'évolution de l'état du système au cours du temps ( entre t k-1 et t k ), une représentation numérique est utilisée. Dans notre cas il s'agit d'un modèle non-linéaire M (SnowModel). Cette représentation se traduit par :

x k = M (x k-1 ) + k-1
(5.1) 94

Théorie de l'assimilation de données

k-1 est l'erreur du modèle au temps t k-1 . Les sources de cette erreur sont nombreuses, mais elles peuvent être rangées en trois catégories :

-Discrétisation du système.

-Simplification de la physique.

-Paramétrage du modèle (les modèles y compris SnowModel ont de nombreux paramètres qui ne sont pas toujours parfaitement connus).

Observations

Les observations sont le deuxième ingrédient nécessaire pour établir un schéma d'assimilation. Elles peuvent provenir de plusieurs sources (par ex. mesures in situ, mesures satellitaires). Elles représentent des quantités physiques qui ne sont pas forcément des sorties du modèle (par exemple dans notre cas nous possédons des observations de surfaces enneigées qui ne sont pas simulées par SnowModel). D'où, la nécessité d'usage d'un opérateur d'observations H permettant de relier les observations aux variables modélisées :

y k = H(x true k ) + o k (5.2) o k représente l'erreur d'observation à l'instant t k . x true k
est l'état exacte du vecteur d'état. Par la suite, pour des raisons de simplification, x, y et feront respectivement référence à x k , y k et k .

Vecteur d'ébauche

La première estimation du vecteur x effectuée à travers un modèle en boucle ouverte M (SnowModel dans notre cas), est le vecteur d'ébauche. Il est noté x b (b désigne background en anglais).

Vecteur d'analyse

Le but de l'assimilation de données est de chercher à représenter l'état d'un système physique x (cf. quantité de SWE stocké dans un manteau neigeux, quantité de précipitations dans une région) de la façon la plus optimale. Le vecteur résultant de cette estimation est nommé analyse et dénoté x a . Notons que chaque vecteur d'analyse à l'instant t k est considéré comme vecteur d'ébauche pour passer de l'instant t k à l'instant t k + 1 . Par la suite, pour des raisons de simplification, x, y et feront respectivement référence à x k , y k et k .

Filtre particulaire

Dans Chap. 5.1.2 et Chap. 5.1.3 étude. Dans cette partie, nous exposons la théorie du filtre particulaire.

L'idée de base du filtre particulaire est de représenter la distribution de probabilité antérieure (prior probability density function, pdf) p(x) au pas de temps discret t k avec un ensemble aléatoire de particules de taille N : x i . La représentation de la distribution de probabilité antérieure p(x) s'écrit sous la forme (van Leeuwen, 2010) :

p(x) = N i=1 1 N δ(x -x i ) (5.3)
δ représente le symbole de Kronecker.

Soit y le vecteur d'observation à l'instant t k . Nous supposons que les erreurs d'observations sont distribuées d'une façon gaussienne. Le théorème de Bayes stipule que :

p(x|y) = p(y|x)p(x) p(y|x)p(x)dx (5.4)
Sachant que :

p(y|x)δ(x -x j )dx = p(y|x j ) (5.5)
La combinaison des Eq. 5.3, Eq. 5.4 et Eq. 5.5 donne :

p(x|y) = 1 N N i=1 w i δ(x -x i ) (5.6)
w i représente le poids de chaque particule comme suit :

w i = p(y|x i ) N j=1 p(y|x j )
(5.7) La densité p(y|x i ) représente la probabilité de densité des observations sachant l'état du modèle x i . Elle est généralement représentée par une fonction gaussienne (van Leeuwen,

p(y|x

i ) = Aexp[- (y -H(x i )) 2 2σ 2 ] (5.8)
σ représente la covariance des erreurs d'observations. A cause du manque de données d'évaluation, elle est généralement prise constante au cours du temps. (y -H(x i ) définit la distance entre les observations et la modélisation, son calcul est détaillé dans le chapitre 5.4.5. A est une constante, sa détermination n'est pas nécessaire, puisque les poids seront normalisés par la suite [START_REF] Blayo | Advanced Data Assimilation for Geosciences : Lecture Notes of the Les Houches School of Physics : Special Issue[END_REF].

Après la phase de calcul du poids, une procédure de sélection est appliquée. Cette dernière consiste à éliminer les particules ayant un poids faible et dupliquer celles qui ont un poids élevé (les détails de cette procédure sont présentés dans le paragraphe 5.4.5.

Résumé de l'article

Dans le chapitre 4 nous avons testé une approche permettant d'estimer le SWE dans le bassin de l'Ourika du Haut Atlas marocain avec un modèle de manteau neigeux distribué. La conclusion de cette étude était de tester l'apport d'un schéma d'assimilation de données afin de réduire les incertitudes du modèle. Dans l'introduction et le chapitre 3, nous avons montré le potentiel de données Sentinel-2 pour mieux contraindre la modélisation du SWE dans le Haut Atlas. Dans l'article ci-dessous, nous présentons donc l'implémentation d'un schéma d'assimilation de données de surfaces enneigées obtenues par la mission satellite Sentinel-2 à très haute résolution spatiale et temporelle dans un modèle de neige alimenté par les forçages météorologiques MERRA-2. Cette assimilation de données est conçue pour corriger les séries de précipitations et de température de l'air fournies par MERRA-2.

Ce schéma d'assimilation de données est basé sur un filtre particulaire. Le modèle est alimenté par les données MERRA-2 comme dans le chapitre précédent. La validation est effectuée en utilisant des observations indépendantes de la hauteur de la neige, du débit et de la surface enneigée (MODIS), mais aussi en utilisant les sorties d'une simulation de référence obtenue en forçant le même modèle avec des données météorologiques in situ enregistrées par les stations automatiques dans le bassin. Les principaux résultats tirés de cette étude sont :

-Les sorties du modèle avec assimilation de données (SWE, P, T) en utilisant MERRA-2 (DA) sont plus proches des sorties de la simulation de référence avec des données in-situ (AWS) que des sorties en utilisant MERRA-2 en boucle ouverte (OL). Une différence notable est la quantité de précipitations en mi-février, où les précipitations estimées par l'assimilation de données sont supérieures à celles estimées à partir des stations météorologiques seulement. Il est difficile de trancher laquelle des simulations est la meilleure avec les données dont nous disposons. -D'après les données MODIS, les trois simulations reproduisent bien l'évolution de la surface enneigée, avec une petite supériorité des résultats de l'assimilation. -La comparaison avec les débits est encourageante malgré l'absence de représentation des transferts hydriques entre la base du manteau neigeux et l'exutoire du bassin versant. Par conséquent il est probable que la calibration d'un modèle hydrologique intermédiaire permettrait d'améliorer cette comparaison.

Assimilation of Sentinel-2 Data into a snowpack

model in the High Atlas mountains

Abstract

The snow melt from the High Atlas is a critical water resource in Morocco. In spite of its importance, monitoring the spatio-temporal evolution of key snow cover properties like the snow water equivalent remains challenging due to the lack of in situ measurements at high elevation. Since 2015, the Sentinel-2 mission provides high spatial resolution images with a 5 day revisit time, which offers new opportunities to characterize snow cover distribution in mountain regions. Here we present a new data assimilation scheme to estimate the state of the snowpack without in situ data. The model was forced using MERRA-2 data and a particle filter was developed to dynamically reduce the biases in temperature and precipitation using Sentinel-2 observations of the snow cover area. The assimilation scheme was implemented using SnowModel, a distributed energy-balance snowpack model and tested in a pilot catchment in the High Atlas. The study period covers 2015-2016 snow season which corresponds to the first operational year of Sentinel-2A, therefore the full

Introduction

The snow melt is an important water resource in many semi-arid and Mediterranean regions [START_REF] Viviroli | Mountains of the world, water towers for humanity : Typology, mapping, and global significance[END_REF]. In the Tensift region near Marrakesh, Morocco, farmers are heavily dependent on melt water from the High Atlas mountains to irrigate the crops during the growing season. Snow melt also provides drinkable water and recharges groundwater [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF]; [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]; Baba et al. (2018a). [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF] estimated that snow melt contribution to the discharge of the five major rivers from the High Atlas mountains ranged from 15% to 45%.

In semi-arid and Mediterranean mountains like the High Atlas, the snow accumulates during winter in the high elevation areas of the catchments, where the terrain complexity enhances the variability of the climatic conditions [START_REF] Fayad | Snow hydrology in Mediterranean mountain regions : A review[END_REF]. Spatial variability in near-surface meteorological variables like precipitation, air temperature, wind speed, solar radiation result in heterogeneous snow accumulation and ablation patterns, which influence the spatio-temporal dynamics of the melting rates and streamflow [START_REF] Marks | A spatially distributed energy balance snowmelt model for application in mountain basins[END_REF]; [START_REF] Fassnacht | Spatio-temporal snowmelt variability across the headwaters of the southern rocky mountains[END_REF].

The main challenge to monitor the snowpack evolution in the High Atlas mountains is the lack of automatic weather stations at high elevation. There is no operational snow course program either and the cost of airborne LIDAR campaigns is prohibitive for water agencies. Satellite remote sensing was used to track the variations of the snow cover area in the High Atlas [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF], but water managers are mainly interested in the snow water equivalent (SWE) and snow melt, which cannot be retrieved by current space-borne sensors [START_REF] Dozier | Estimating the spatial distribution of snow water equivalent in the world's mountains[END_REF].

Global-scale climate reanalyzes like MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2 Gelaro et al. (2017b)) provide atmospheric data that are continuous in space and time. These data can be used instead of weather station data to simulate the snow water equivalent in any mountain region, provided that a spatial downscaling is applied to account for first-order effects of surface topography on key variables like air temperature and precipitation Durand et al. (2008a); [START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF]; Baba et al. (2018a). In the High Atlas, MERRA-2 data enabled to simulate reasonably well the evolution of the snow cover over [START_REF] Whiteman | Mountain meteorology : fundamentals and applications[END_REF][START_REF] Africanews | Morocco experiences heavy snowfall, social media reacts[END_REF][START_REF] Baba | Contribution of high resolution remote sensing data to the modeling of the snow cover the in Atlas Mountains[END_REF]. However, some large biases were identified in the precipitation forcing. Indeed, although current climate reanalyzes integrate many satellite-based and conventional weather observations, their output is still subject to significant biases, especially in mountainous regions and in regions with few meteorological stations like Africa [START_REF] Gelaro | The modern-era retrospective analysis for research and applications, version 2 (MERRA-2)[END_REF].

Complementary remote sensing observations can be used to avoid the spreading of forcing errors in the model through data assimilation approaches. In particular, the snow cover area is easily retrieved from multi-spectral optical imagery [START_REF] Dietz | Remote sensing of snow-a review of available methods[END_REF]; [START_REF] Dumont | Optical remote sensing of snow cover[END_REF]. Given the wide availability of snow cover products from MODIS or Landsat, many studies have focused on the assimilation of remotely sensed snow cover area (or snow cover fraction) into a snowpack model to dynamically adjust the model parameters or to correct errors in model forcing data [START_REF] Andreadis | Assimilating remotely sensed snow observations into a macroscale hydrology model[END_REF]; [START_REF] Clark | Hydrological data assimilation with the ensemble Kalman filter : Use of streamflow observations to update states in a distributed hydrological model[END_REF]; Durand et al. (2008a); [START_REF] Zaitchik | Forward-looking assimilation of MODIS-derived snow-covered area into a land surface model[END_REF][START_REF] De Lannoy | Multiscale assimilation of advanced microwave scanning radiometer-EOS snow water equivalent and moderate resolution imaging spectroradiometer snow cover fraction observations in northern Colorado[END_REF]; [START_REF] Thirel | Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter[END_REF]; [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF]; [START_REF] Stigter | Assimilation of snow cover and snow depth into a snow model to estimate snow water equivalent and snowmelt runoff in a Himalayan catchment[END_REF]; [START_REF] Toure | Assimilation of MODIS snow cover fraction observations into the NASA catchment land surface model[END_REF]. In particular, assimilation of snow cover area data improved the simulation of SWE in mountain regions with partial or transient snow cover, whereas the improvement was marginal in regions with deep snowpacks [START_REF] Toure | Assimilation of MODIS snow cover fraction observations into the NASA catchment land surface model[END_REF]; [START_REF] Zaitchik | Forward-looking assimilation of MODIS-derived snow-covered area into a land surface model[END_REF]. The assimilation of high resolution snow cover area images from Landsat into a distributed energy balance model enabled to significantly reduce the bias in SWE simulation in the California Sierra Nevada and the subtropical Andes [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF]; [START_REF] Cortés | Impacts of El Niño and La Niña on interannual snow accumulation in the andes : Results from a high-resolution 31 year reanalysis[END_REF]. Landsat data were preferred to MODIS data in both cases, based on the premise that only Landsat data can accurately capture the spatial variability of the snow cover in semi-arid and Mediterranean mountains [START_REF] Margulis | A particle batch smoother approach to snow water equivalent estimation[END_REF]. However, a limitation of the Landsat mission for snow cover applications is its long revisit time of 16 days. Given the reported cloud probability in Mediterranean mountain regions [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]; [START_REF] Gascoin | A snow cover climatology for the Pyrenees from MODIS snow products[END_REF], the probability to have a clear observation with Landsat-8 can be lower than once per month per pixel. Since 2017, the Sentinel-2 mission provides high resolution multi-spectral observations of the global land surface with a 5-day revisit time [START_REF] Gascon | Copernicus Sentinel-2A calibration and products validation status[END_REF]. Hence, Sentinel-2 data have the potential to further improve data assimilation schemes for snow modelling. In addition, as part of the Copernicus program of the European Commission, Sentinel-2 data are freely available and should be provided continuously for the next decade at least [START_REF] Drusch | Sentinel-2 : ESA's optical highresolution mission for GMES operational services[END_REF], which are two important assets in the perspective of an operational snow monitoring system in a developing country like Morocco.

The objective of this study is to implement an assimilation scheme for Sentinel-2 data to simulate the spatial distribution of the snow water equivalent and snow melt in the High Atlas without ground data. We introduce a new assimilation scheme based on Snow-Model, a spatially-distributed energy-balance snowpack model Liston and Elder (2006a), which has been previously evaluated in the High Atlas [START_REF] Baba | Contribution of high resolution remote sensing data to the modeling of the snow cover the in Atlas Mountains[END_REF]. The data assimilation scheme is based on the particle filter, as this approach is better adapted to non-linear physically-based model like SnowModel [START_REF] Moradkhani | Uncertainty assessment of hydrologic model states and parameters : Sequential data assimilation using the particle filter[END_REF]; [START_REF] Leisenring | Snow water equivalent prediction using Bayesian data assimilation methods[END_REF]; van [START_REF] Van Leeuwen | Particle filters for applications in geosciences[END_REF]. The model takes MERRA-2 data as meteorological forcing and Sentinel-2 snow cover area products as observational data. The other input data are a digital elevation model and a land cover map. The scheme is demonstrated in the case of the Rheraya catchment over the 2016-2017 snow season. The Rheraya catchment is a pilot catchment in the High Atlas near Marrakech where meteorological data are available for the evaluation of the data assimilation [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. We compare the output of the data assimilation scheme to the open loop simulation (i.e., without assimilation of Sentinel-2 data), and to a reference simulation which was forced by in situ meteorological data instead of MERRA-2.

To our knowledge, this study is the first to use Sentinel-2 data to estimate SWE in mountain regions. It is also the first attempt to assimilate remote sensing data to characterize snow distribution in the High Atlas.

Study area and data

Study area

The study area is the Rheraya catchment in the High Atlas Mountains of Morocco (Fig. 5.4). The catchment elevation ranges from 1000 m (river gauge station of Tahanaout) to 4167 m (Mount Toubkal, the highest peak in the North Africa). The Rheraya river discharge at Tahanaout is strongly influenced by the snow melt contribution in spring (Fig. 5.5). The discharge during the melt season between March and May represents about 50 % of the annual discharge. The catchment is sparsely vegetated, except in the valleys bottoms, where snowfalls rarely occur.

We selected the Rheraya catchment because it is the only catchment in the High Atlas with high elevation automatic weather stations [START_REF] Jarlan | Remote sensing of water resources in semi-arid mediterranean areas : the joint international laboratory TREMA[END_REF]. The Rheraya catchment was also the study site of previous snow hydrology studies in Morocco [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF][START_REF] Marchane | Climate change impacts on surface water resources in the Rheraya catchment (High Atlas, Morocco)[END_REF][START_REF] Hajhouji | Evaluation d'un modèle pluie-débit sur le bassin versant de la Rheraya dans le Haut Atlas marocain[END_REF].

Data Topography and land cover

The 90 m Shuttle Radar Topography Mission (SRTM) digital elevation model was used after resampling it to 200 m spatial resolution by using cubic interpolation. The choice of the spatial resolution is based on a previous study in the Rheraya catchment [START_REF] Baba | Effect of digital elevation model resolution on the simulation of the snow cover evolution in the High Atlas[END_REF], where we analyzed the trade-off between spatial resolution and the computation time. This study suggested that a distributed snow model with a grid coarser than 250 m is not able to capture the high spatial variability in the snow cover patterns due to the importance of the incoming solar radiation in the snowpack energy budget.

The land cover use data set is derived from European Space Agency Climate Change I (ESA-CCI). The majority of the area is bare soil.

Snow cover area

Sentinel-2 We queried all available Sentinel-2 snow covered area (SCA) products between January and July 2016 in the L2B-Snow collection of the Theia Land data center [START_REF] Gascoin | Theia snow collection[END_REF]. All snow products were derived from Sentinel-2A acquisitions since Sentinel-2B was only launched in 2017. The Rheraya catchment intersects tiles 29RNQ and 29RPQ, but we used only the products from tile 29RPQ since the tile 29RNQ only covers the low elevation area of the catchment where the snow is rarely observed. The downloaded snow products were clipped over the Rheraya catchment in their native projection (WGS-84 UTM 29N). To reduce the input/output operation time in the assimilation phase, we removed the products with a cloud cover higher than 70%. We obtained 13 SCA maps of the Rheraya catchment with a snow cover fraction ranging from 82% to 0% and a cloud fraction ranging from 0% to 66% (Table 5.2).

These SCA products were generated from Sentinel-2 images based on the MAJA/LIS processor. The MAJA processor computes the surface reflectance and the cloud and cloud shadow mask [START_REF] Lonjou | Maccs-atcor joint algorithm (MAJA)[END_REF]). The output of MAJA is a level-2A product, that is read by the LIS processor to determine the snow cover area at 20 m resolution [START_REF] Simon | Algorithm theoretical basis documentation for an operational snow cover product from Sentinel-2 and Landsat-8 data (let-it-snow[END_REF]. The snow detection is performed using the "flat surface reflectances", i.e. surface reflectances that were corrected to remove the first order effect of the topography. The snow classification in LIS uses the Normalized Difference Snow Index [START_REF] Dozier | Spectral signature of alpine snow cover from the Landsat Thematic Mapper[END_REF] :

NDSI = ρ green -ρ SWIR ρ green + ρ SWIR , ( 5.9)
where ρ green (resp. ρ SWIR ) is the flat surface reflectance in the green channel (resp. SWIR at 1.6 µm). The NDSI encapsulates the fact that almost only snow surfaces are very bright in the visible but dark in the shortwave infrared. The output of LIS is a single band raster resolution having four possible values : snow, no-snow, cloud and no-data. For further details, the reader can refer to the Algorithm Theoretical Basis Documentations of MAJA [START_REF] Hagolle | MAJA algorithm theoretical basis document[END_REF] and LIS [START_REF] Simon | Algorithm theoretical basis documentation for an operational snow cover product from Sentinel-2 and Landsat-8 data (let-it-snow[END_REF].

MODIS We used the daily MODIS snow products (MOD10A1 collection 6) to validate the different simulations. These products have a spatial resolution close to 500 m and a [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]Baba et al., 2018a). First, the NDSI in each snow pixel was converted to SCF using the equation that was used in the collection 5 [START_REF] Hall | Algorithm theoretical basis document (ATBD) for the modis snow and sea ice-mapping algorithms[END_REF] SCF = -0.01 + 1.45 × NDSI (5.10)

Then, the data gaps due to cloud obstruction were interpolated with the algorithm of [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF], which was specifically developed and validated in the High Atlas mountains. The products were finally clipped over the Rheraya catchment.

MERRA-2 data

The meteorological forcing were taken from NASA's second Modern-Era Retrospective analysis for Research and Applications (MERRA-2). MERRA-2 reanalysis data are widely used, including studies dealing with the assimilation of snow cover data [START_REF] Baldo | Towards large-scale implementation of a high resolution snow reanalysis over midlatitude montane ranges[END_REF]; [START_REF] Cortés | Impacts of El Niño and La Niña on interannual snow accumulation in the andes : Results from a high-resolution 31 year reanalysis[END_REF]; [START_REF] Toure | Assimilation of MODIS snow cover fraction observations into the NASA catchment land surface model[END_REF]. MERRA-2 data have a spatial resolution of 0.5 • × 0.625 • and a temporal resolution of 1 hour. MERRA-2 includes the assimilation of satellite observations that were not available to MERRA and an improved representation of key variables for simulating snowpack evolution [START_REF] Gelaro | The modern-era retrospective analysis for research and applications, version 2 (MERRA-2)[END_REF]. We extracted the precipitation, 2-m air temperature, 2-m specific humidity and 2-m wind speed in the four nearest grid cell of the Rheraya catchment. The post-processing steps to make MERRA-2 compatible with SnowModel are described in Baba et al. (2018a).

In situ observations

The hourly data from three automatic weather stations (AWS) were used as input to SnowModel to generate synthetic observations and as validation data (snow height at Oukaimeden). Table 5.3 summarizes the AWS characteristics and available data. We also used a continuous record of the snow height at Oukaimeden (3230 m.a.s.l.) measured by an acoustic snow gauge.

We also used the observed river discharge at the gauging station of Tahanaout (31.29 o N,-7.96 o E, 1048 m) provided by the Tensift river basin water agency (ABHT). A river stage device measures the height of the water three times per day. The gauging is performed once per month at best for logistical reasons. These measurements are known to have large errors given the difficulty to establish a robust stage-discharge relation in the Rheraya river, which has an unstable river bed.

Methods

Validation strategy

We focused on the 2015-2016 snow season. We defined three simulations :

-the open loop simulation (OL) was obtained by running SnowModel with downscaled MERRA-2 data as input. The output of this simulation is called the prior. -the data assimilation simulation (DA) was obtained in the same configuration as the open loop simulation but downscaled MERRA-2 forcings were perturbed to assimilate the Sentinel-2 snow cover maps through a particle filter. The output of this simulation is called the posterior. -the synthetic data simulation (referred to as AWS) was obtained by running Snow-Model with in situ meteorological observations from the AWS. The output of this simulation was considered as an independent dataset to evaluate the effect of the assimilation.

The main input and output of the simulations is presented in Fig. 5.6.

The performance of the DA was qualified by computing the RMSE (Eq. 5.11) and the Pearson correlation coefficient between the output of the DA scheme and the synthetic observations of air temperature, precipitation, SWE and SCA.

RMSE = 1 N N i=1 (X i -X truth i ) 2 (5.11) R = E[XY ]E[X]E[Y ] E[X 2 ] -E[X] 2 E[Y 2 ] -E[Y ] 2
(5.12)

where X i is the ith element of the vector X, and X truth i is the ith element of the synthetic observation vector X truth . E is the expectation. Y is the observation vector.

The output of the OL, DA and AWS were also compared to in situ observations of the snow height at Oukaimeden, the snow cover fraction of the catchment area from MODIS data and the river discharge at Tahanaout. Regarding the river discharge comparison, the challenge is that SnowModel does not simulate the water transfer from the snow melt to the streamflow. The runoff is the sum of the rain and melt water released from the snowpack. Given the complexity of the intermediate hydrological processes in the Rheraya catchment [START_REF] Chaponnière | Understanding hydrological processes with scarce data in a mountain environment[END_REF][START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF], we compared the observed discharge with the simulated runoff using 3-day averages. The simulated discharge was the description of MicroMet submodel only.

In MicroMet, a DEM is used to correct the effect of elevation on air temperature, humidity and precipitation. In the case of the air temperature, the following equation is used :

T x = T stn -τ (Z x -Z stn ) (5.13)
where T stn is the air temperature at the station, Z stn the station elevation, T x ( o C) the air temperature in the target grid cell, Z x the elevation of the grid cell. τ is the temperature lapse rate ( o C.km -1 ), which depends on the month of the year. As a first step, MicroMet adjusts the station data to a common elevation Z = 0 using Eq. 5.13. Then, the temperature is interpolated using a Gaussian distance-dependent weighting function to the model grid [START_REF] Barnes | A technique for maximizing details in numerical weather map analysis[END_REF]. Finally, Eq. 5.13 is used again to adjust the gridded temperatures values back to the actual elevation taken from the DEM.

The same method is used for the relative humidity and precipitation. The relative humidity at the station is first converted to dew point temperature. The same function as for the air temperature is used to account for elevation (Eq 5.13), but with different lapse rate values. In the case of the precipitation, a different function of elevation is used : .14) where P stn is the precipitation at the station, P x (mm) the precipitation in the target grid cell, Z x the elevation of the grid cell. χ is the precipitation correction factor (km -1 ), which also depends on the month of the year.

P x = P stn × [1 + χ(Z x -Z stn )] [1 -χ(Z x -Z stn )] ( 5
The shortwave and longwave radiations are derived from the interpolated temperature and humidity on each grid cell accounting for the effect of the slope and orientation of the terrain [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF].

As noted above, MicroMet was originally designed to spatially distribute meteorological variables from weather station measurements. However it can also be used to downscale large scale meteorological data like MERRA-2 data. In this case, the MERRA-2 data are simply considered as the meteorological measurements of a virtual weather station that is located in the center of the MERRA-2 grid cell. The elevation of the virtual station is calculated from the surface geopotential height of the MERRA-2 cell. This approach was successfully applied in the High-Atlas Baba et al. (2018a) and in the Andes [START_REF] Mernild | The Andes Cordillera. part i : snow distribution, properties, and trends (1979-2014)[END_REF].

Preliminary sensitivity tests revealed that the parameterization of the precipitation solid fraction caused significant biases in the simulated snow cover area after the snowfalls. Excessive model biases should be corrected as much as possible to make the assimilation scheme efficient. By default, MicroMet uses a temperature threshold of 2 o C to distinguish between solid and liquid precipitation [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. We found that theses biases could be reduced by enhancing the default parameterization in SnowModel to account for the influence of the relative humidity on the precipitation phase using the empirical relation of [START_REF] Froidurot | Sensitivity of precipitation phase over the Swiss Alps to different meteorological variables[END_REF] :

F s = 1 1 + e γRH+βT +α (5.15)
where F s is the solid fraction in the total precipitation, RH is the relative humidity in % and T the 2 m air temperature during the same time step. The parameters α, γ, β were optimized on a different snow season. We chose the 2008-2009 snow season for which we dispose of an accurate time series of high resolution SCA maps from Formosat-2 data [START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]. A set of 200 (α, γ, β) triplets were generated in the intervals [22.00, 26.00], [-2.80, -2.40], [-0.25, -0.20] by Monte-Carlo sampling and used as parameters of the same SnowModel run as in [START_REF] Baba | Effect of digital elevation model resolution on the simulation of the snow cover evolution in the High Atlas[END_REF]. These intervals were chosen following the recommendations of [START_REF] Froidurot | Sensitivity of precipitation phase over the Swiss Alps to different meteorological variables[END_REF]. The optimization was done using the same model configuration as presented in [START_REF] Baba | Contribution of high resolution remote sensing data to the modeling of the snow cover the in Atlas Mountains[END_REF]. We used the Heidle Skill Score (HSS) [START_REF] Notarnicola | Snow cover maps from MODIS images at 250 m resolution, part 2 : Validation[END_REF] between the simulated snow cover area and the Formosat-2 SCA maps over the entire study period to identify the optimal parameters triplet (α = 25, β = -2.5,γ = -0.2).

SWE-SCA conversion

We define the SCA as a binary variable indicating the presence (SCA=1) or absence (SCA=0) of snow in a grid cell, and the SCF as a continuous variable indicating the snow covered fraction of a grid cell. The SCA or the SCF are not standard outputs of SnowModel hence must be derived from the simulated SWE or snow depth. This transformation is done through a "snow depletion curve". There are several snow depletion curve formulations in the literature [START_REF] Anderson | Twentieth-century drought in the conterminous United States[END_REF][START_REF] Liston | Representing subgrid snow cover heterogeneities in regional and global models[END_REF], and their parameters are poorly known especially in the context of the High Atlas. In addition snow depletion curves were mostly developed for low to mid-resolution imagery like MODIS. Therefore, in this study we chose to assimilate the SCA rather than the SCF to reduce the importance of the snow depletion curve equation and parameters.

First, we compute the SCF at the grid cell level using the snow depletion curve of [START_REF] Zaitchik | Forward-looking assimilation of MODIS-derived snow-covered area into a land surface model[END_REF], also used by [START_REF] Thirel | Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter[END_REF] :

SCF = min(1.0, 1 -[exp( -τ SWE SWE SCF=1 ) - SWE SWE SCF=1 exp(-τ )]) (5.16)
where SWE is the snow water equivalent in the grid cell, τ is the snow distribution shape parameter relating the total amount of snow in the grid cell to the percent snow cover within the grid cell (by default it is set to 4.0 [START_REF] Thirel | Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter[END_REF]). SWE max defines the minimum SWE to reach a full snow coverage within the grid cell, whose value depends on the land cover. In this study we used the value for bare soil (SWE max 13 mm). The grid cell was considered as snow covered if the SCF is greater than a fixed value SCF 0 , i.e. : SCA = u[SCF -SCF 0 ]

(5.17)

where u is the discrete form of the Heaviside function. SCF 0 was set to 0.25.

Data assimilation algorithm

In this study, the data assimilation scheme is used to reduce the biases in air temperature and precipitation. Let x the state vector containing a discrete approximation of the SWE, air temperature, precipitation. The basis of the particle filter is to seek the best estimation of the prognostic variables, by generating an ensemble of N random draws (particles) x i where i ∈ [1 : N ]. These particles are propagated through SnowModel whenever an observation (here a Sentinel-2 snow product) is available. At each observation date, a specific normalized weight is assigned to every particle [START_REF] Van Leeuwen | Nonlinear data assimilation in geosciences : an extremely efficient particle filter[END_REF][START_REF] Van Leeuwen | Particle filters for applications in geosciences[END_REF] :

w i,j = 1 2σ √ 2π e -( i,j ) 2 2σ 2
( 5.18) where w i,j is the weight of the ith particle at the jth observation date, σ is the standard deviation of the observation error. We assume that σ is constant throughout the study period. Here sigma was set to 0.15 based on the comparison with very high resolution images in the Alps and a visual inspection of the snow products in our study area. i,j is often defined as the difference between the spatial average of the observed and modeled SCF [START_REF] Leisenring | Snow water equivalent prediction using Bayesian data assimilation methods[END_REF]. Here, we aim to take into account the differences between the model and the observations at the grid cell level (i.e. on a pixel basis) when defining the quality of each particle. Hence we chose to derive i,j from the Heidle skill score (HSS), a statistical index based on the confusion matrix, which is often used to evaluate the skills of snow remote sensing products at the pixel level (Notarnicola et al., 2013a) :

HSS = 2(T P × T N -F P × F N ) (T P + F P )Ö(F P + T N ) + (T P + F N )(F N + T N ) (5.19)
where T P , T N , F P , F N are the number of true positive, true negative, false positive and false negative pixels, respectively. The perfect simulation has an HSS equal to 1 while the worst has an HSS equal to 0. Hence, i,j is defined as :

i,j = 1 -HSS i,j .

(5.20)

Figure 5.7 -Principle of the stochastic universal sampling method in the case of 7 particles [START_REF] Pencheva | Modelling of a stochastic universal sampling selection operator in genetic algorithms using generalized nets[END_REF]. The brown numbers identify each ith particle and the segment size is proportional to its normalized weight W i .

Subsequently, each weight is updated using

w * i,j = n obs j=1
w i,j (5.21) where w * i,j is the updated value of the weight, n obs is the number of assimilated observations until this time step (including the jth). Thereafter the weights are normalized using :

W i = w * i,j N i=1 w i (5.22)
where W i is the updated normalized weight. After computing the weights, the particle filter removes particles with low weight and duplicates particles with high weight to estimate the posterior distribution of the state vector. This step is called resampling. Many resampling algorithms exist (see a review by [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF]). The most used algorithm in hydrological studies is the stochastic universal sampling (SUS) method, due to the ease of its implementation, its low computational time and also its satisfying results [START_REF] Noh | Applying sequential Monte Carlo methods into a distributed hydrologic model : lagged particle filtering approach with regularization[END_REF][START_REF] Pencheva | Modelling of a stochastic universal sampling selection operator in genetic algorithms using generalized nets[END_REF].

First, a vector containing the cumulative sum of W i is computed. Then, N pointers spaced by 1/N are used to select the particles in this space (Fig. 5.7). This enables to preferentially select a particle with a high weight given that this particle will span a larger segment of this space, hence particles with the highest weight will be duplicated. The selected particles are thereafter perturbed as explained in Sect. 5.4.5.

An issue arises if all particles have equivalent weights, e.g. if the domain is fully snow covered or snow-free. In this case, every particle is equally likely to be selected for the next assimilation window, which will prevent the creation of new particles since the number of particles must remain constant. Therefore, we modified the SUS to select, at each assimilation date, only 50% of the particles for the next assimilation window, i.e. we use only N SU S = N/2 pointers instead of N . Thus, each selected particle is duplicated at least Fig. 5.8 shows the difference between the standard SUS resampling and the enhanced resampling in such a case of equivalent particle weights. In the standard algorithm, only one particle can be perturbed in the next step, whereas in the enhanced algorithm, three particles can be perturbed. This will increase the probability to obtain a better estimation of the state vector for the next assimilation step.

To summarize, we present the different steps of the computational implementation of the DA scheme as follow, and also illustrate them in Fig. 5.9.

1. Sample N ens particles x i by perturbing MERRA-2 precipitation and temperature.

2. Integrate all particles in SnowModel forward time (from t to t + 1).

3. Calculate the weights according to Eq. 5.22. 6. Choose the best particle. The model output (SWE) that correspond to the most likely SW E true state.

Ensemble of meteorological forcings

At each assimilation step, the MERRA-2 air temperature and precipitation are perturbed to generate an ensemble of meteorological forcing spanning the next assimilation window. It is important to emphasize that the perturbation is done before the downscaling by MicroMet, i.e. from the coarse scale MERRA-2 data, which are considered as virtual station measurements (Sect. 5.4.5).

For the air temperature, the perturbation is done by adding a white noise with zero mean and standard deviation of 2 o C [START_REF] Clark | Hydrological data assimilation with the ensemble Kalman filter : Use of streamflow observations to update states in a distributed hydrological model[END_REF][START_REF] Leisenring | Snow water equivalent prediction using Bayesian data assimilation methods[END_REF]. The precipitation data were perturbed by a multiplicative correction factor randomly sampled between 0.75 and 1.5 [START_REF] Charrois | On the assimilation of optical reflectances and snow depth observations into a detailed snowpack model[END_REF].

Implementation

The DA workflow uses code in Matlab, Fortran and shell scripts. The scheme was implemented in the CNES high performing computing cluster (HPCC). The main script uses the Portable Batch System (PBS Pro) job scheduler to execute the different steps. As a first step, the main HPCC PBS script creates N ens perturbed meteorological forcing files (ASCII) from the initial meteorological file. Thereafter, N ens simulations are performed with SnowModel in parallel (Figure 5.10). This allows to get N ens of SWE simulations (GDAT binary files). For each particle the SWE is converted to SCA, then its weight is computed based on the observed snow cover area (GeoTiff). Finally, all particles enter

Results

Comparison to in situ snow height

The different simulations captured the overall snow height evolution at Oukaimeden (Fig. 5.11), however the DA simulation outperformed the OL and AWS simulation in terms of correlation and RMSE (Tab. 5.4).

There remain large discrepancies from the end of February to the mid-March even in the DA simulation. Over the same period there is no assimilation data. In the three cases the melting rates were too high especially in March. By contrast, the DA simulation captured well the snow height peaks on 16 February and 26 March. Both peaks occurred between two successive Sentinel-2A acquisitions (07 and 17 February, 18 and 28 March, Tab. 5.2). Comparison to MODIS Fig. 5.12 shows the evolution of the snow covered fraction of the Rheraya catchment over the study period from the three simulations and from MODIS observations. The simulated snow covered fraction of the catchment was obtained by taking the average of the SCF in every grid cell, which was obtained by applying Eq 5.16. The three simulations reproduce well the MODIS observations (Tab 5.4). The DA cSCF is larger than the OL cSCF, meaning that the correction of the precipitation led to increase the snowfall. The RMSE and correlation both indicate that the DA has improved the realism of the simulation with respect to the OL simulation (Tab. 5.4). The DA simulation is close to the AWS forcing but it is also slighlty better in terms of RMSE and correlation.

Comparison to discharge observations

Despite the crude approach to compute the discharge from the model output, the three simulations provide a realistic representation of the discharge at the 3-day time step. The DA simulation better captures the spring melt flood at the end of March (Fig. 5.13). In addition, the AWS simulation has a large error in mid-April, which is not present in the OL and DA simulations. This error is likely due to wrong precipitation data in the AWS record.

Comparison to the AWS simulation

As explained in Sect 5.4.5, we use the AWS simulation to assess the effect of the data assimilation on the SWE, air temperature and precipitation. Figure 5.14 shows the evolution of the mean SWE over the Rheraya catchment from the three simulations. The difference between the DA SWE and OL SWE is highest between mid-February to mid-March. The DA has led to a significant increase in the simulated SWE with respect to the OL simulation. As a result, the DA SWE is closer to the AWS SWE during the same period. As shown in Figure 5.16, the higher SWE in the DA run is the result of higher precipitations, in particular the DA has increased the precipitation rate of the two major events on 14 February and 27 March. However, despite these differences in the accumulation period, the DA and OL simulations produced a similar SWE in the end of the melt season. This can be explained by the differences in evolution of the air temperature over the melt season (Fig. 5.15). The mean air temperature in the DA simulation is higher than the one of the OL simulation between February to May, which increased the melting rate in the DA simulation.

On 05 January, we note that the DA has reduced the OL precipitation (Fig 5.16), which is consistent with the AWS simulation (Fig 5.14).

Discussion

The results show that the data assimilation scheme has reduced the bias and increased the correlation with two independent observation datasets of snow height and snow cover area. However, these observations provide only a partial evaluation of the simulations : (i) the snow height was measured only at the Oukaimeden weather station, therefore it is not necessarily representative of the snow height evolution over the entire catchment, even at similar elevations ; (ii) a better agreement with the MODIS snow cover area does not necessarily means that other snowpack properties like the SWE were better simulated. However, the combination of these observations suggest that the assimilation had a positive influence on the model outputs.

The comparison with the synthetic data (AWS simulation) showed that MERRA-2 provided a good estimate of the meteorological forcing in agreement with a previous study in the neighbouring catchment of the Ourika catchment (Baba et al., 2018a). This means that there was a relatively low bias in the prior precipitation and air temperature, which could have caused otherwise a divergence in the posterior SWE. Yet, the assimilation generally tended to increase the precipitation and SWE (Fig. 5.16, Fig. 5.17). However, it decreased the SWE in the highest area of the catchment near the Toubkal peak (Fig. 5.17. This is consistent with the fact that MicroMet uses a crude representation of the orographic enhancement of the precipitation (Sect. 5.4.5). The equation 5.14 can lead to unrealistic of the precipitation at high elevations if the reference stations used in the interpolation are located at much lower elevations. Here the elevation of MERRA-2 virtual stations ranges between 836 m and 1892 m, whereas the elevations in southern part of the catchment exceed 3000 m.a.s.l. During the melt season, the data assimilation also increased the air temperature (Fig. 5.15) and as a result the air temperature is more consistent with the AWS simulation. Therefore it means that the posterior air temperature Figure 5.17 -Map of the differences between the mean SWE from simulations DA and OL. The mean SWE was computed over the entire study period. is more consistent with in situ observations of the air temperature.

Figure 5.18 shows the evolution of the HSS with respect to the Sentinel-2 data over the simulation period. As expected, the HSS of the DA simulation is higher than the OL HSS since the HSS was used in the function to weight the particles during the assimilation. The median of the HSS values is 0.72 for OL, 0.81 for DA and 0.71 for AWS. The improvement is mostly significant on 17 February and 07 April. The assimilation of the 17 February Sentinel-2 image in particular has caused an increase in the SWE which persisted until mid-March (Fig. 5.14). This image was captured only 3 days after the precipitation event of the 14 February, the largest precipitation event of the study period (Fig. 5.16). As a result this observation was important in the general performance of the assimilation scheme.

The approach presented in this study is based on the downscaling of large scale MERRA-2 data with MicroMet. Each MERRA-2 grid cell is considered as a virtual station. In the proposed data assimilation scheme, the meteorological forcing is perturbed at each virtual station before downscaling to the resolution of the snowpack model. This enables to restrain the number of particles (and thus the computation cost), in comparison with a data assimilation scheme which perturbs the forcing at the model grid cell level. Another advantage of this approach is that it guarantees a spatially smooth distribution of the posterior SWE, because MicroMet distributes the meteorological forcing using an optimal Gaussian distance-dependent weighting function [START_REF] Barnes | A technique for maximizing details in numerical weather map analysis[END_REF]; Liston and Elder, ). However, a drawback of this method is that localized weather phenomenons like convective precipitation events might be not be correctly corrected. Another important limitation of our approach is related to the generation of the perturbed precipitation data.

Here we used a stochastic multiplicative factor, which means that if the MERRA-2 precipitation is zero, the perturbed precipitation data is unchanged. Hence, the success of the particle filter is strongly dependent upon the quality of the prior precipitation from MERRA-2. In this study, we found that MERRA-2 provided a good first guess for the snow season 2015-2016. A previous study in the High Atlas showed that the general quality of MERRA-2 data is higher over the most recent years, probably due to the increasing amount of analyzed weather satellite data (Baba et al., 2018a).

Conclusions

We presented a new data assimilation scheme to integrate Sentinel-2 snow cover area in the simulation of the SWE in the High Atlas. The scheme is a particle filter with a modified stochastic universal sampling. The snowpack model is the distributed snow evolution model SnowModel. The method is based only on remote sensing and global climate data hence it does not require ground measurements. The results showed that the assimilation has improved the simulation of the snow cover, although this study covered only the 2015-2016 snow season before the launch of Sentinel-2B when the 5-day revisit time of the Sentinel-2 mission was not yet achieved. Since late 2017, both Sentinel-2 satellites are acquiring data over all land masses at full systematic 5-day revisit. This should increase the value of Sentinel-2 observations for snow studies in particular in areas of transient snow cover like semi-arid and Mediterranean climate mountain areas. Indeed, our study showed that a single image could have a significant positive impact on the posterior SWE, e.g. if it is acquired shortly after a snowfall event.

Here we focused on the assimilation of the snow cover area because it is easily retrieved from widely available optical remote sensing observations. However, assimilating the snow cover area only can cause lead to ambiguous posterior estimations. For instance, if the prior underestimates the snow cover duration at high elevation with respect to satellite observations, this can be corrected by decreasing the air temperature (to reduce melting rate) but also by increasing the precipitation (to increase snow accumulation). This equifinality issue can be tackled by incorporating additional observations. In particular, we argue that even sporadic SWE or snow height measurements (e.g. snow courses) would drastically reduce the uncertainties in the posterior SWE. If such measurements are not performed by operational weather and water agencies, an option is to use crowd-sourced snow data [START_REF] Giuliani | Using crowdsourced web content for informing water systems operations in snow-dominated catchments[END_REF]. Our data assimilation scheme can be adapted to integrate this kind of observations because it is based on a distributed snow modelling. Additional remote sensing products like wet snow areas from Sentinel-1 could also be incorporated. Finally, additional meteorological variables could be updated using the same assimilation scheme, e.g. the shortwave radiation, which is an important driver of the snowpack energy-balance in the High Atlas [START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF]. In this perspective, further model developments are also required to enhance the representation of the radiation processes in SnowModel, including the effect of Saharan dust deposition on snow albedo.

The continuous improvement of global-scale reanalyses and the increasing availability of remote sensing data should lead to a better characterization of the SWE in the High Atlas even without ground measurements. Ultimately, this should contribute to a better knowledge of the water resources for the local communities. d snow is omnipresent between December and March. Chapitre 6

Conclusion

Synthèse

L'eau issue de la fonte du manteau neigeux dans le Haut Atlas marocain est une ressource vitale pour l'agriculture dans les zones en aval. Pourtant, peu de stations météorologiques sont installées dans les bassins versants de haute altitude, ce qui limite notre capacité à suivre l'état du manteau neigeux au cours de la saison.

L'objectif principal de notre travail consistait à proposer des méthodes adaptées au contexte marocain pour suivre cette ressource naturelle. Pour parvenir à cet objectif nous avons été amenés à effectuer trois études principales. Les résultats de ces études ont été présentés sous la forme de trois articles publiés dans différentes revues scientifiques.

Dans un premier temps (Chap. 3), nous avons cherché la résolution spatiale qui permet de simuler l'évolution du manteau neigeux de façon distribuée sans que l'effet de la topographie soit perdu d'une part, et que le temps de calcul du modèle avec cette résolution soit raisonnable d'autre part. Pour arriver à trouver cette résolution optimale, nous avons généré un MNT à très haute résolution spatiale (4 m) à partir de trois couples d'images stéréoscopiques Pléiades. Ce MNT a été ré-échantillonné à différentes résolutions allant de 8 m à 1000 m pour examiner l'effet de la résolution du modèle sur la simulation du manteau neigeux au cours de la saison 2008-2009 dans le bassin de la Rheraya avec le modèle de manteau neigeux SnowModel. Nous avons évalué la qualité de chaque simulation à l'aide d'une série d'images Formosat-2 à 8 m de résolution. Il ressort de cette étude que dans le bassin de la Rheraya, un MNT de 250 m présente le meilleur compromis recherché (qualité de simulation/temps d'exécution). Ce résultat a été corroboré à l'aide d'une analyse complémentaire sur le semivariogramme des pentes, ce qui permet d'envisager une méthode indépendante pour le transposer dans d'autres régions.

Dans un deuxième temps (Chap. 4), nous avons évalué l'apport des données météorologiques MERRA-2 pour l'estimation de l'état du manteau neigeux dans le Haut Atlas sur une plus longue période (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). La modélisation du manteau neigeux a été effectuée avec SnowModel à une résolution de 250 m comme recommandé précédemment. A la lumière des résultats obtenus, nous avons conclu qu'il est possible de produire une estimation réaliste de l'évolution du manteau neigeux par cette méthode malgré la résolution grossière des données MERRA-2. Cela est possible grâce aux algorithmes de spatialisation des forçages météorologiques adaptés aux régions de montagne qui sont implémentés dans la routine MicroMet Liston and Elder (2006a). Cependant, la qualité des simulations est variable selon les années. En effet, les performances sont médiocres pour certaines années (surtout en début de période). Ceci est principalement dû à des erreurs dans les forçages météorologiques (particulièrement les précipitations) qui se sont propagées au cours de la saison. Ainsi, nous avons recommandé de corriger les biais potentiels dans les données météorologiques MERRA-2 en utilisant un schéma d'assimilation de données basé sur les surfaces enneigées obtenues par satellite comme proposé dans la littérature [START_REF] Cortés | Snow process estimation over the extratropical andes using a data assimilation framework integrating MERRA data and Landsat imagery[END_REF].

Dans un troisième temps (Chap. 5), nous avons conçu, implémenté, et évalué un schéma d'assimilation de données de type filtre particulaire. Ce schéma a été alimenté par des données de surfaces enneigées à haute résolution dérivées de Sentinel-2A pour la saison 2015-2016 dans le bassin de la Rheraya. Les résultats de cette étude suggèrent que la combinaison des données MERRA-2, des observations satellitaires de surfaces enneigées et d'un modèle distribué de neige, à travers un schéma d'assimilation basé sur le filtre particulaire, permet de bien estimer l'état du manteau neigeux. La méthode fonctionne sans mesures in situ et s'appuie sur des données gratuites. Elle est donc transposable sur n'importe quel autre bassin versant du Haut Atlas voire d'autres massifs montagneux.

Discussion

Nous notons que pour la saison où nous avons testé le schéma d'assimilation, seules les données Sentinel-2A étaient disponibles. Désormais, nous pouvons aussi compter sur les données de Sentinel-2B qui est opérationnel depuis 2017. Ainsi la résolution temporelle sera réduite de 10 jours à 5 jours (Fig. 6.1), ce qui augmentera le nombre d'images de la surface enneigée et donc devrait contribuer à améliorer les performances du schéma d'assimilation.

L'évaluation des simulations dans les différents chapitres se base principalement sur des observations de surfaces enneigées, de la hauteur de la neige à Oukaimeden (3230 m.a.s.l) et des débits à l'exutoire des deux bassins versants. Un autre type de comparaison a été utilisé dans le chapitre 5, qui est basé sur la génération d'observations synthétiques à partir des observations météorologiques issues des stations. Cette simulation a été considérée comme la plus réaliste pour représenter le SWE dans le bassin de la Rheraya. Or, aucune validation de l'équivalent en eau ou de la hauteur de la neige n'a été effectuée de façon distribuée à partir de données d'observation. Ainsi, les conclusions de cette étude sont à prendre avec prudence. Pour bien évaluer les erreurs du modèle et de l'assimilation, il paraît nécessaire de confronter les sorties du modèle avec des données de références distribuées. Un autre problème lié à notre méthode d'assimilation concerne l'équifinalité1 des estimations postérieures de température et des précipitations. Il est possible que deux chroniques différentes de précipitation et de température permettent de bien simuler l'état observé de la surface enneigée. Seules des mesures complémentaires météorologiques ou du SWE pourraient permettre de lever ce type d'ambiguïté quand elle se présente.

Enfin, il est important de rappeler que les différentes études effectuées durant cette thèse se sont concentrées sur le versant ouest du Haut Atlas qui est exposé à l'influence climatique de l'océan Atlantique, par opposition au versant qui est davantage exposé à l'influence du climat désertique du Sahara. Il est probable que les conditions plus arides sur le versant oriental favorisent un taux de sublimation plus élevé. [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF] ont trouvé que ce taux variait de 7% à 20% (Fig. 6.2) entre 2003 et 2010 à Oukaimeden dans le bassin de la Rheraya. Par ailleurs, [START_REF] Schulz | Snowmelt and sublimation : field experiments and modelling in the High Atlas Mountains of Morocco[END_REF] ont estimé que 44% de la neige disparaît par sublimation dans la montagne de M'Goun, c'est-à-dire dans la partie Est du Haut Atlas, à cause des longues périodes sèches et froides sur ce site. Sur cette base, il serait intéressant de prolonger ce travail à partir des données MERRA-2 pour évaluer la variabilité des taux de sublimation. Pour cela un effort supplémentaire de validation des données de vitesse du vent devrait être mené avec des anémomètres sur le terrain.

Recommandations

Compte tenu de l'importance de la neige comme ressource en eau, et des fortes incertitudes qui persistent dans la caractérisation du manteau neigeux malgré les avancées réalisées au cours de cette thèse, nous recommandons de créer un réseau météorologique [START_REF] Boudhar | Energy fluxes and melt rate of a seasonal snow cover in the Moroccan High Atlas[END_REF]. pérenne dans le Haut Atlas marocain. Certes, les résultats de la télédétection et de la modélisation sont encourageants, mais l'intégration de données météorologiques et de hauteurs de neige in situ pour la validation ou la réduction des erreurs du modèle est nécessaire pour établir un système efficace d'aide à la décision.

Nous recommandons aussi de tester le schéma d'assimilation dans une autre zone où une base de données de validation plus complète est disponible. Le bassin versant de la Tuolumne en Californie pourrait être un bon site d'étude puisque la NASA fournit des mesures LIDAR de la hauteur de la neige tous les quinze jours environ sur plusieurs saisons [START_REF] Painter | The airborne snow observatory : Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo[END_REF]. Ainsi, ces nouveaux jeux de données nous permettraient de valider précisément et spatialement les sorties du schéma d'assimilation. The simulated SWE at every Formosat-2 acquisition date was converted to binary snow maps using a SWE 0 threshold. Then, the simulated snow maps were resampled to the same grid as Formosat-2. This enables to compute a confusion matrix between both images, from which we can derive the statistics described in Sect.

Figure 1 . 1 -

 11 Figure 1.1 -Enneigement et localisation des montagnes du Haut Atlas dans la région méditerranéenne.

1. 1 .

 1 Figure 1.2 -La fonte de la neige est primordiale pour l'alimentation en eau en aval. Photo prise à Tachedirt au Haut Atlas (Décembre 2016).

Figure 1 . 4 -

 14 Figure 1.4 -Schéma décrivant les principaux éléments de l'équation du bilan de masse et d'énergie.

Figure 1 . 5 -

 15 Figure 1.5 -Gradients mensuels optimisés de température de l'air (G) pour la saison de 2008-2009 pour le bassin de la Rheraya. R 2 représente son coefficient de détermination par rapport aux observations de la température issues par des stations météorologiques.

Figure 1 . 6 -

 16 Figure 1.6 -Photographie du processus de carottage.

Chapitre 1 -Figure 1 . 7 -

 117 Figure 1.7 -Structure d'un modèle de neige

Figure 1 .

 1 Figure 1.10 -Comparaison visuelle des images MODIS et Sentinel-2 au Haut-Atlas Marocain. Les zones enneigées sont déterminées avec meilleure précision avec Sentinel-2.

Figure 1 .

 1 Figure 1.11 -Réflectance spectrale de la neige et des nuages en fonction de la longueur d'onde dans l'intervalle 0-3µm. La barre rouge correspond à la bande B11 infrarouge de Sentinel-2.

Figure 1 .

 1 Figure 1.12 -Exemple d'une image Formosat-2 du bassin de la Rheraya. La distinction neige/nuage peut se faire à l'oeil.

Figure 1 .

 1 Figure 1.13 -Photographie prise aux montagnes du Haut Atlas, montrant la forte variabilité spatiale de la distribution de la couverture neigeuse. Décembre 2016. Méthode Variables Résolution spatiale

Figure 1 .

 1 Figure 1.14 -Illustration mettant en exergue les apports et les limitations de la télédétection et de la modélisation menant à la proposition d'implémentation d'un schéma d'assimilation de données. Les phrases en vert représentent les apports, et celles en rouge représentent les limitations.

Figure 2 . 1 -

 21 Figure 2.1 -Bassin versant de Tensift et ses sous-bassins atlasiques (B.V.). Sa position géographique dans la carte du Maroc (point rouge) et la localisation des bassins versants. L'étoile rouge représente la ville de Marrakech. Les exutoires sont présentés par des ronds bleues et leurs noms sont en bleu.

Figure 2 . 2 -

 22 Figure 2.2 -Précipitations saisonnières moyennes (1963-2003) pour les différentes stations météorologiques.

2. 1 .Figure 2 . 3 -

 123 Figure 2.3 -Débits moyens (m 3 /s) pour chaque saison. Aghbalou : exutoire de l'Ourika, Tahanaout : exutoire de la Rheraya, Iguir N'Kouris exutoire de N'Fis, Taferiat exutoire de Zat et Sidi Rahal exutoire de R'dat.

Figure 2 . 4 -

 24 Figure 2.4 -Précipitations et débits moyens saisonniers observés à Aghbalou.

Figure 2 . 5 -

 25 Figure 2.5 -Carte de répartition des précipitations moyennes annuelles dans le bassin versant du Tensift[START_REF] Boudhar | Télédétection du manteau neigeux et modélisation de la contribution des eaux de fonte des neiges aux débits des Oueds du Haut Atlas de Marrakech[END_REF].

Figure 2 .

 2 Figure 2.7 -Processus des interactions des composantes de SnowModel utilisées dans cette thèse

Figure 2 . 8 -

 28 Figure 2.8 -Simplification du processus de spatialisation de la température par Micro-Met. De (1) à (2) : Utilisation du facteur de correction de la température (temperature lapse rate) pour passer les températures observées des altitudes des stations, aux températures à la surface de référence (Z = 0 m). De (2) à (3) Spatialisation des températures à la surface de référence à l'aide du schéma de Barnes. De (3) à (4) Deuxième utilisation du facteur de correction de la température (temperature lapse rate) (Eq. 2.4)

Figure 2 .

 2 Figure 2.10 -Chaîne de traitement des images stéréo-Pléiades afin de générer un MNT à très haute résolution spatiale. 33

Figure 2 .

 2 Figure 2.11 -Les différents produits de télédétection utilisés durant cette thèse. Formosat-2, Sentinel-2 et MODIS ont permis d'extraire les surfaces enneigées durant les saisons mentionnées. Tandis que Pléiades a permis de générer un modèle numérique terrain à très haute résolution spatiale.

Figure 2 .

 2 Figure 2.12 -Photographie de la station météorologique d'Oukaimeden.

Figure 2 .

 2 Figure 2.13 -Évolution de la SCF en fonction du SWE avec les deux approches : a) Seuil fixe, b) SDC de Zaitchik and Rodell (2009) (Eq. 2.10).

Figure 2 .

 2 Figure 2.14 -Principe de détermination de la matrice de confusion et du calcul de ses éléments.

Figure 3 . 1 -

 31 Figure 3.1 -Photograph taken near Oukaimeden in the Rheraya catchment on 14 December 2016. (31.194 o N, -7.584 o S). It illustrates the effect of slope and aspect on the snow cover area at decametric scales that is often observed in the High Atlas.

Figure 3 . 2 -

 32 Figure 3.2 -Map of the study area. The map on the left shows the digital elevation model derived from Pléiades stereo-images over the Rheraya catchment, the extent of the Formosat-2 acquisitions, the location of the AWS.

Figure 3 . 3 -

 33 Figure 3.3 -Characteristics of the Pléiades images. For each site, base-to-height ratio (B/H) the ratio of the distance between two successive positions of the satellite to its height above ground, is an indicator of the sensitivity to topography.

Figure 3 . 4 -

 34 Figure 3.4 -Comparison of the control points (red dots) and the ortho-image.

Figure 3 . 5 -

 35 Figure 3.5 -Times series of Formosat-2 images used in this study. The SCA derived from the images is represented in green polygons and superposed to natural color composites of the images.

Figure 3 . 6 -

 36 Figure 3.6 -Flow chart of the method used to compute the performance of the model at different resolution based on the Formosat-2 snow maps.The simulated SWE at every Formosat-2 acquisition date was converted to binary snow maps using a SWE 0 threshold. Then, the simulated snow maps were resampled to the same grid as Formosat-2. This enables to compute a confusion matrix between both images, from which we can derive the statistics described in Sect.3.2.4.

Figure 3 .

 3 Figure 3.7 -Time series of the snow cover area at each resolution with the default precipitation lapse rates (top) and the optimized precipitation lapse rates (bottom). The Formosat-2 observations are shown in black circles. The snow cover area is given in fraction of the Rheraya catchment area.

Figure 3 .

 3 Figure 3.8 -Comparison of the HSS between the simulated and the Formosat-2 snow maps at every resolution. The boxplots show the distribution of HSS for every SWE 0 values and every date (see Sect. 3.2.4 for details on the method).

Figure 3 .

 3 Figure 3.9 -Comparison of the observed number of snow days (Formosat-2) and the simulated number of snow days for each model resolution (SM) . The number of snow days were computed using the 15 dates of Formosat-2 acquisitions. All maps are shown with the same colorbar in units of days.

Figure 3 .

 3 Figure 3.10 -Simulated and measured snow depth at Oukaimeden (From November 2008 to June 2009) for every model resolution.

Figure 3 .

 3 Figure 3.11 -Spatial averages of a) the model grid elevation, b) the model grid slope, c) the total snow precipitation, d) the mean SWE, e) the mean incoming shortwave radiation, f) the mean daily melt, for each model resolution. These values were extracted from the validation area, i.e. the intersubsection of the Formosat-2 extent and the catchment area.

Figure 3 .

 3 Figure 3.12 -Simulated SWE on 01 April 2009. The black line indicates the extent of the snow cover from the Formosat-2 image acquired on the same day. The red line indicates the model domain boundary.

  3.14). The semivariogram computed in the interval 0-600 m further shows that the break is located at 200 m (inset in left panel of Fig.3.14). In addition, the slope of the semivariogram in the double logarithmic plot has a break near 200 m, which indicates a change in the fractal dimension near this distance.

Figure 3 .

 3 Figure 3.13 -Mean simulated shortwave radiation in the same part of the catchment as Fig. 3.12.

Figure 3 .

 3 Figure 3.14 -Empirical semivariograms (γ) of the slope from the 30 m ASTER GDEM of the study area (validation domain). The semivariograms were computed using 20 bins of distance (h in meters) and a maximum distance of 2000 m. The inset in the left panel shows the semivariogram with a maximum distance of 600 m. The right panel shows the double logarithmic representation of semivariogram, and was computed using logspaced distance values instead of linearly spaced distance values. Every semivariogram was computed using a sub-sample of 30000 randomly drawn points.

Figure 4 . 1 -

 41 Figure 4.1 -Location of the catchments of the five main tributaries of the Tensift river in the High Atlas. SD : Oukaimeden station snow depth measurements.

Figure 4 . 2 -

 42 Figure 4.2 -Mean monthly discharge of the Ourika river at Aghbalou (1970-2003 (Agence Japonaise de coopération internationale, 2008)).

Figure 4 . 3 -

 43 Figure 4.3 -Ourika land cover map.

Figure 4 . 4 -

 44 Figure 4.4 -Summary of the methodology.

  4.1) (Oukaimeden-CAF 2612 m.a.s.l, Agouns 2381 m.a.s.l and Amenzal 2374 m.a.s.l). These data are available from 2000 to 2009.

Figure 4 . 5 -

 45 Figure 4.5 -Time series of the simulated (blue) and observed (red) daily snow depth (SD) at Oukaimeden.

Figure 4 . 6 -

 46 Figure 4.6 -Time series of the simulated (blue) and observed (red) snow cover fraction of the Ourika catchment.

Figure 4 . 7 -Fig. 4 .

 474 Figure 4.7 -Evolution of monthly precipitation between 1 st December to 1 st April from rain gauge measurements (Obs) and MERRA-2/MicroMet (Mod).

Figure 4 .

 4 Figure 4.8 -Time series of the monthly simulated runoff in m 3 /s (blue) and observed discharge in m 3 /s (red) at Aghbalou between January and May from 2000 to 2017.

Figure 4 .

 4 Figure 4.9 -Time series of 10-day step simulated runoff in m 3 /s (blue) and observed discharge in m 3 /s (red) at Aghbalou between January and May from 2008 to 2017.

Figure 4 .

 4 Figure 4.10 -Evolution of the simulated SWE in the Ourika catchment from September 2000 to April 2018. Each snow season from 01 September 2000 to 31 May 2017 was plotted at the daily time step as an individual gray line. The blue line shows the mean of these values and the gray envelope indicates their first and third quantiles. The red line corresponds to the hydrological year 2017-2018. 2014-2015 year is the green one,2006-2007 is the pink and 2004-2005 is the black. Doted lines represent the interval 80 mm -100 mm.

Figure 4 .

 4 Figure 4.11 -From top to bottom : annual values of the maximum simulated SWE, the date of maximum SWE and D.date represents the snow disappearance date (first day after 01-March that has SWE=0).

Figure 4 .

 4 Figure 4.12 -Maps of the simulated 01 April SWE for every year of the study period.

  Ce chapitre se présente sous forme d'un article publié à Remote Sensing : Baba, M. W.,Gascoin, S., and Hanich, L. Assimilation of Sentinel-2 Data into a snowpack model in the High Atlas mountains(V.F https: // www. mdpi. com/ 2072-4292/ 10/ 12/ 1982)

Figure 5 . 1 -

 51 Figure 5.1 -Le filtre de Kalman compare les prédictions des lois de Newton avec les mesures embarquées (cercles rouges) afin de générer une meilleure estimation de la position réelle de l'engin spatial (ligne bleue avec tirets bleus).

Figure 5 . 2 -

 52 Figure 5.2 -Exemple d'analyse de Cressman. L'état de l'ébauche (modèle sans correction) est représenté par la courbe grise pointillée, les observations par les points verts et l'état analysé (après correction) par la courbe rouge.

5. 3 .

 3 Figure 5.3 -Le schéma de gauche représente la modélisation en boucle ouverte, et celui de droite décrit une simulation d'ensemble de particules. Cette dernière consiste à générer des perturbations des conditions initiales afin de réduire les erreurs de modélisation ou de forçages.

Figure 5 . 4 -

 54 Figure 5.4 -Elevations of the Rheraya catchment, and the location of discharge measurement (sky blue), in situ snow height (HS) and automatic weather stations (AWS). The dashed points refers to the intersection between Sentinel-2 images borders and the Rheraya.

Figure 5 . 5 -

 55 Figure 5.5 -Mean monthly discharge (m 3 /s) of the Rheraya River at Tahanaout (river discharge data from the Agence du Bassin Hydraulique du Tensift (1974-2003).

Figure 5 .

 5 Figure 5.8 -Schema describing the weight calculation and resampling. (1) : with initial SUS resampling. (2) with enhanced SUS method.

Figure 5 .

 5 Figure 5.9 -Description of the DA scheme used in this study

Figure 5 .

 5 Figure 5.10 -Description of the IT implementation of the PF in this study.

Figure 5 .

 5 Figure 5.11 -Comparison of the simulated and observed in situ snow height in Oukaimeden at the daily time step. The blue diamonds in the upper part of the figure indicate the date of the assimilated Sentinel-2 observations.

Figure 5 .

 5 Figure 5.12 -Comparison of the simulated and observed snow covered fraction of the Rheraya catchment (cSCF) at the daily time step. The blue diamonds in the upper part of the figure indicate the date of the assimilated Sentinel-2 observations.

Figure 5 .

 5 Figure 5.13 -Comparison of the simulated and observed discharge at Tahanaout (see Sect. 5.4.5)

Figure 5 .

 5 Figure 5.14 -Comparison of prior, posterior and AWS modeled SWE.

Figure 5 .

 5 Figure 5.15 -Evolution of the mean air temperature over the Rheraya catchment during i) all season ii) melt season.

Figure 5 .

 5 Figure 5.16 -Comparison between the OL and DA mean daily precipitation over the Rheraya catchment.

Figure 5 .

 5 Figure 5.18 -Temporal evolution of HSS (compared to Sentinel-2 SCA maps) for OL, DA and AWS simulations.

  2006b

Figure 6 . 1 -

 61 Figure 6.1 -Évolution de la moyenne du SWE entre 2000-2018 sur le bassin de l'Ourika avec une approximation du temps de passage des satellites Sentinel-2A et Sentinel-2B combinés durant une année hydrologique (les tirets verticaux).

Figure 6 . 2 -

 62 Figure 6.2 -Variations du taux de sublimation (Échelle entre 0% et 20%) par année à Oukaimeden (Boudhar et al., 2016).

3. 2

 2 Map of the study area. The map on the left shows the digital elevation model derived from Pléiades stereo-images over the Rheraya catchment, the extent of the Formosat-2 acquisitions, the location of the AWS. . . . . 46 3.3 Characteristics of the Pléiades images. For each site, base-to-height ratio (B/H) the ratio of the distance between two successive positions of the satellite to its height above ground, is an indicator of the sensitivity to topography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Comparison of the control points (red dots) and the ortho-image. . . . . . 49 3.5 Times series of Formosat-2 images used in this study. The SCA derived from the images is represented in green polygons and superposed to natural color composites of the images. . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6 Flow chart of the method used to compute the performance of the model at different resolution based on the Formosat-2 snow maps.

  3.2.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.7 Time series of the snow cover area at each resolution with the default precipitation lapse rates (top) and the optimized precipitation lapse rates (bottom). The Formosat-2 observations are shown in black circles. The snow cover area is given in fraction of the Rheraya catchment area. . . . . 57 3.8 Comparison of the HSS between the simulated and the Formosat-2 snow maps at every resolution. The boxplots show the distribution of HSS for every SWE 0 values and every date (see Sect. 3.2.4 for details on the method). 58 3.9 Comparison of the observed number of snow days (Formosat-2) and the simulated number of snow days for each model resolution (SM) . The number of snow days were computed using the 15 dates of Formosat-2 acquisitions. All maps are shown with the same colorbar in units of days. . . . . . . . . 3.10 Simulated and measured snow depth at Oukaimeden (From November 2008 to June 2009) for every model resolution. . . . . . . . . . . . . . . . . . . . 3.11 Spatial averages of a) the model grid elevation, b) the model grid slope, c) the total snow precipitation, d) the mean SWE, e) the mean incoming shortwave radiation, f) the mean daily melt, for each model resolution. These values were extracted from the validation area, i.e. the intersubsection of the Formosat-2 extent and the catchment area. . . . . . . . . . . . 3.12 Simulated SWE on 01 April 2009. The black line indicates the extent of the snow cover from the Formosat-2 image acquired on the same day. The red line indicates the model domain boundary. . . . . . . . . . . . . . . . . 3.13 Mean simulated shortwave radiation in the same part of the catchment as Fig. 3.12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Empirical semivariograms (γ) of the slope from the 30 m ASTER GDEM of the study area (validation domain). The semivariograms were computed using 20 bins of distance (h in meters) and a maximum distance of 2000 m. The inset in the left panel shows the semivariogram with a maximum distance of 600 m. The right panel shows the double logarithmic representation of semivariogram, and was computed using log-spaced distance values instead of linearly spaced distance values. Every semivariogram was computed using a sub-sample of 30000 randomly drawn points. . . . . . . . . . . . . 4.4 R 2 between the simulated and observed snow cover area (SCA) and discharge (Q) in the Ourika catchment for every hydrological year of the study period. The R 2 was computed at the daily timestep from 01 September to 31 May for SCA and at the monthly timestep from January to May for Q (NA : no measurements available for this year). . . . . . . . . . . . . 5.1 Notations utilisées dans ce chapitre. . . . . . . . . . . . . . . . . . . . . . . 5.2 Acquisition dates of Sentinel-2 snow products and fraction of the Rheraya catchment covered by snow and clouds. . . . . . . . . . . . . . . . . . . . . 5.3 Description of the three automatic weather stations (T : temperature, P : precipitation, RH : relative humidity, HS : height of snow) . . . . . . . . . 5.4 RMSE and correlation between modeled and observed (i) daily snow height (HS) and (ii) snow cover fraction over the catchment area (cSCF) from 01 September to 01 May. HS was measured at Oukaimeden and cSCF was derived from MODIS data. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 1 -Caractéristiques des satellites/capteurs permettant de suivre la couverture neigeuse gratuitement en temps quasi-réel à l'échelle globale.

	Chapitre 1 -Introduction
	Satellite	Capteur Résolution spatiale (m) Temps de revisite Fauchée
	Landsat 8	OLI	15-30	16 j	185 km
	Terra-Aqua MODIS 500	1 j	2330 km
	Sentinel-2	MSI	10-20	5 j	290 km

  ρ vis et ρ mir représentent la réflectance dans le visible et la réflectance dans le moyen infrarouge respectivement. Un pixel est considéré comme enneigé lorsque le NDSI est supérieur à un seuil (généralement entre 0.2 et 0.4).

			.
	Il est défini par la formule suivante :
	NDSI =	ρ vis -ρ mir ρ vis + ρ mir	(1.7)

-Classification supervisée : cette méthode consiste à lever manuellement un jeu de données d'apprentissage qui permet d'entraîner un algorithme de classification. Cette méthode a été utilisée par

[START_REF] Marchane | Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range[END_REF]

pour classifier les zones enneigées dans une série d'images Formosat-2 à 8 m de résolution, car ce capteur étant dépourvu d'un canal dans le moyen infrarouge, la méthode non-supervisée du NDSI ne peut pas être utilisée (Fig.

1

.12). Même une méthode supervisée ne permet pas toujours de distinguer la neige des nuages en l'absence de moyen infrarouge. Un examen visuel complémentaire est recommandé pour éviter les confusions.

Table 1 .

 1

2 -Récapitulatif des différents d'approches du suivi de l'enneigement l'estimation du SWE dans un bassin versant. Dans cette section, nous présenterons l'apport de chacune d'elles pour le suivi du manteau neigeux au Haut Atlas du Maroc, ainsi que ses limitations, afin de trancher pour la méthode qui sera la plus optimale à notre contexte. Les montagnes du Haut Atlas ont un climat semi-aride, et sont caractérisées par une forte variabilité spatio-temporelle des précipitations et de la température. De ce fait, la distribution de la neige est très hétérogène (Fig.

1

.13) et sera difficilement capturée par un modèle de neige semi-distribué ou conceptuel. Les modèles distribués pourront à leur tour être d'un grand intérêt mais ils nécessiteront des forçages météorologiques fiables afin de tourner sans trop être entachés par les erreurs. Les rares mesures in situ existantes ne permettent pas d'échantillonner correctement la variabilité spatio-temporelle du climat et de la neige à l'échelle du Haut Atlas (Chap.1.2.2). Le tableau. 1.4 et la figure. 1.14 récapitulent les apports et limitations de chaque approche.

stratégie envisagée dans cette thèse est donc d'implémenter un schéma d'assimilation de surfaces enneigées dans un modèle de neige.

  1.4. Problématique de suivre l'évolution du manteaux neigeux. En effet, comme le résume le tableau. 1.4 ces deux approches sont complémentaires, puisque la télédétection offre une estimation de la couverture neigeuse à large échelle, et la modélisation permet d'extraire les variables clés pour l'hydrologie. Mais la modélisation est entachée de fortes incertitudes liées aux erreurs de forçages et à la complexité des processus physiques imparfaitement représentés par les modèles même les plus sophistiqués. Pour éviter la propagation des erreurs de forçage dans le modèle, la Ce schéma doit permettre de corriger les erreurs de forçages météorologiques afin d'avoir une meilleure estimation de l'état du manteau neigeux et donc une meilleure connaissance du SWE et de la fonte principalement.

Pour réaliser ce travail, notre étude sera décomposée en trois parties complémentaires : -Résolution spatiale : Le schéma d'assimilation proposé est basé sur la réalisation d'ensembles de simulations distribuées, ce qui implique de limiter tant que possible le temps d'exécution de chaque simulation. Nous avons donc, dans un premier temps, cherché la résolution spatiale optimale permettant de faire tourner ce modèle distribué sans trop dégrader les simulations spatialisées (Chapitre. 3). -Modélisation en boucle ouverte : Nous avons appliqué ce même modèle de neige distribué en l'alimentant avec des forçages des réanalyses MERRA-2 afin de déterminer sa capacité à reproduire l'état du manteau neigeux pour une période s'étalant de 2000 à 2018 (Chapitre. 4). Les résultats de ce chapitre plaident en faveur de l'implémentation d'un schéma d'assimilation de données basé sur MERRA-2. -Assimilation de données : Le chapitre 5 consiste à concevoir, implémenter et valider un schéma d'assimilation de surface enneigées dans un modèle de neige pour estimer la distribution spatiale de l'équivalent en eau.

  'une superficie de 20456 km 2 , la plaine du Haouz est drainée par l'oued Tensift et ses affluents venus du Haut Atlas (Fig.2.1) et est caractérisée principalement par les cultures de céréales, d'oliviers et d'abricotiers. Le sud de cette plaine est marqué par la présence de la chaîne des montagnes du Haut Atlas dont les sommets dépassent fréquemment 3000 m et qui culmine à 4167 m au djebel Toubkal. Les principaux sous-bassins versants de cette chaîne montagneuse sont : N'Fis, Rheraya, Ourika, Zat et R'Dat. Les caractéristiques de ces sous-bassins à l'exutoire de la station de jaugeage sont présentées dans la table2.1.1.

	Chapitre 2
	Outils, méthodes et zones d'études
	2.1 Zones d'études
	2.1.1 Région

D

Table 2 .

 2

	2.5). Les précipitations se

1 -Caractéristiques des sous-bassins versants atlasiques du bassin versant Tensift (bassins jaugés). Altitude max. et altitude moy. correspondent respectivement à l'altitude maximale et moyenne.

Table 2 .

 2

	Nom de la station Altitude (m) Données météorologiques Bassin
	Aremd	1965	T, RH, P, RS	Rheraya
	Imskerbour	1404	T, RH, P	Rheraya
	Oukaimeden	3230	T, RH, P, HS	Rheraya/Ourika
	Tamatert	1924	P	Rheraya
	Oukaimeden-CAF 2612	P	Ourika
	Neltner	3207	T,RH,P	Rheraya
	Agouns	2387	P	Ourika
	Amenzal	2374	P	Ourika
	2 -Caractéristiques des différentes stations de mesures de la zone d'études. T :
	température, RH : humidité relative (relative humidity en anglais). P : précipitation, HS :
	hauteur de la neige(snow height). RS : rayonnement solaire.

Trois pluviomètres sont installés dans le bassin versant d'Ourika (le plus haut est à 2387 m d'altitude). Dans le bassin de la Rheraya, plusieurs stations météorologiques et pluviomètres sont installés (dont deux stations sont localisées au delà de 3000 m). Ces stations mesurent différentes variables météorologiques. La table. 2.1.3 résume les

Figure 2.6 -Localisation des données hydro-météorologiques utilisées au cours de cette étude. Oukaimeden HS est la station de mesure de la hauteur de la neige (snow height).

Saison 2008-2009

	Pour cette saison, nous possédons une série de quinze cartes d'enneigement issues
	d'images du satellite Formosat-2 à 8 m de résolution. Elles couvrent toute la saison de fonte
	(entre le 8 février et le 10 juin). Pour faciliter l'extraction des surfaces enneigées et éviter
	les sources d'erreurs liées à la présence des nuages, seules les images à faible couverture
	nuageuses ont été sélectionnées. L'extraction des surfaces enneigées a été effectuée en
	utilisant une classification supervisée par

Effet de la résolution spatiale sur la modélisation du manteau neigeux dans le Haut Atlas

				Chapitre 3
	Ce chapitre se présente sous forme d'un article publié à Water Resources Research :
				).
	Dans le cas inverse, la formule se transforme en :
	lim d→∞	HSS =	2TP 2TP + FP + FN	(2.13)
	Dans cette étude, nous considérons une classe comme prédominante lorsque sa surface
	est supérieure à 4 fois celle de l'autre classe.

Table 3 .

 3 2 -Precipitation correction factors (χ in km -1 ) and temperature lapse rates (τ in o C.m -1 ) in the study area.Lapse rate : Jan Feb Mar Apr May Jun Jul Aug Sep Oct

	Nov Dec

Table 3 .

 3

	3 -Snowmodel parameters
	Snow Model parameters	Value Unit
	Curvature length scale	500	m
	Slope weight	0.58	-
	Curvature weight	0.42	-
	Melting snow cover albedo 0.6	-
	Dry snow albedo	0.8	-
	Initial snow density	0.3	-

Table 3 .

 3 4 -RMSE of the simulated snow cover area (SCA, in fraction of the catchment area) with respect to Formosat-2 observations and RMSE of the simulated snow depth (SD) with respect to Oukaimeden in situ measurements.

	Spatial resolution 8 m 30 m 90 m 250 m 500 m 1000 m
	SCA (%)	6.5 3.9	5.7	12.1	16.7	16.9
	SD (cm)	30.2 30.6 30.0 30.1	47.4	64.3
	with the model grid of 250 m (300 model realizations were done).

Table 4 .

 4 1 -Characteristics of the five main tributaries of the Tensift river in the High Atlas. Z 0 : elevation of the river gage station, z : mean elevation of the catchment, Q : mean annual discharge, Q S : mean annual specific discharge, F S : mean annual fraction of snow melt in the discharge. The discharge are reported for the period 1970-2003 (Agence Japonaise de coopération internationale, 2008). The fraction of snow melt was computed over[START_REF] Hiemstra | Snow redistribution by wind and interactions with vegetation at upper treeline in the Medicine Bow Mountains, Wyoming, USA[END_REF] -2005[START_REF] Boudhar | Evaluation of the snowmelt runoff model in the Moroccan High Atlas Mountains using two snow-cover estimates[END_REF].

	mm) F S

Table 4 .

 4 2 -MERRA-2 variables used in this study.

	Variable name	Description	Unit
	PHIS	Geopotential height m 2 /s 2
	LON	Longitude	o
	LAT	Latitude	o
	T2M	Temperature at 2m	K
	QV2M	2m specific humidity kg.kg -1
	U2M	2m eastward wind	m.s -1
	V2M	2m northward wind m.s -1
	PRECTOTCORR Precipitation rate	kg.m -2 .s -1

Table 4 .

 4  bias and correlation between simulated and observed snow depth at

	Oukaimeden (m)
	Year	04-05 05-06 08-09 10-11 13-14 14-15
	RMSE 0.21	0.19	0.22	0.15	0.33	0.17
	Bias	0.06	-0.07 0.08	0.03	-0.11 0.01
	R 2	0.37	0.88	0.82	0.49	0.66	0.92

Table 5 .

 5 nous avons présenté un aperçu général des différentes méthodes d'assimilation et nous avons justifié le choix du filtre particulaire pour notre 1 -Notations utilisées dans ce chapitre.

	Symbole Définition	Remarques
	x	Vecteur d'état
	y	Vecteur d'observation
	x b	Vecteur d'ébauche	Première estimation du vecteur d'état
	x a	Vecteur d'analyse	Le plus proche de la réalité
	M	Modèle
	H	Opérateur d'observations
	t	Temps
	t k	Temps discret
	o	Erreur d'observation

Table 5 .

 5 2 -Acquisition dates of Sentinel-2 snow products and fraction of the Rheraya catchment covered by snow and clouds.

	Dates	08-01-16 18-01-16 07-02-16 17-02-16 18-03-16 28-03-16 07-04-16
	SCA(%)	1.19	1.94	1.54	79.93	19.72	35.83	18.11
	Cloud(%)	63.33	13.69	0.75	1.50	7.02	3.33	2.45
	Dates	27-04-16 07-05-16 27-05-16 06-06-16 26-06-16 06-07-16
	SCA(%)	4.16	9.09	0.11	0.23	0.05	0.0
	Cloud(%)	1.87	4.33	41.28	0.00	69.13	0.0

Table 5 .

 5 3 -Description of the three automatic weather stations (T : temperature, P : precipitation, RH : relative humidity, HS : height of snow)

	Stations	Coordinate (WGS 84)	Elevation (m) Available data
	Neltner	(31.063 o N, -7.938 o E) 3207	P,T,RH
	Tachedirt	(31.158 o N, -7.849 o E) 2393	P,T,RH
	Imskerbour	(31.205 o N, -7.938 o E) 1404	P,T,RH
	Oukaimeden (31.180 o N, -7.865 o E) 3230	HS
	daily revisit time. The MODIS products were post-processed to generate a gap-free time
	series of snow cover fraction images (SCF)

Table 5 .

 5 4 -RMSE and correlation between modeled and observed (i) daily snow height (HS) and (ii) snow cover fraction over the catchment area (cSCF) from 01 September to 01 May. HS was measured at Oukaimeden and cSCF was derived from MODIS data.

		OL	DA AWS
	RMSE (HS)	13.65 9.08 14.36
	R (HS)	0.75 0.82 0.64
	RMSE (cSCF) 8.25 6.40 8.70
	R (cSCF)	0.87 0.92 0.86

Table des figures

 des 1.1 Enneigement et localisation des montagnes du Haut Atlas dans la région méditerranéenne. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 La fonte de la neige est primordiale pour l'alimentation en eau en aval. Photo prise à Tachedirt au Haut Atlas (Décembre 2016). . . . . . . . . . . 3 1.3 Différentes formes de précipitations selon la température de l'air. Au niveau des hautes altitudes où la température est négative, les précipitations sont sous forme solide (cf. neige). Dans les zones où la température est élevée, les précipitations sont sous forme liquide. Alors que dans les zones où la température est proche de zéro, les deux différentes formes de précipitations sont possibles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Schéma décrivant les principaux éléments de l'équation du bilan de masse et d'énergie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Gradients mensuels optimisés de température de l'air (G) pour la saison de 2008-2009 pour le bassin de la Rheraya. R 2 représente son coefficient de détermination par rapport aux observations de la température issues par des stations météorologiques. . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.6 Photographie du processus de carottage. . . . . . . . . . . . . . . . . . . . 11 1.7 Structure d'un modèle de neige . . . . . . . . . . . . . . . . . . . . . . . . 12 1.8 Exemple des trois approches de modélisation appliquées à une zone du bassin de la Rheraya. Les points représentent la modélisation ponctuelle. Les fenêtres colorées représentent la modélisation semi-distribuée. Les petits carreaux représentent les mailles de la modélisation distribuée. . . . . . . . 13 1.9 Illustration de l'apport de la télédétection spatiale pour le suivi de l'enneigement dans les zones montagneuses. . . . . . . . . . . . . . . . . . . . . . 14 1.10 Comparaison visuelle des images MODIS et Sentinel-2 au Haut-Atlas Marocain. Les zones enneigées sont déterminées avec meilleure précision avec Sentinel-2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.11 Réflectance spectrale de la neige et des nuages en fonction de la longueur d'onde dans l'intervalle 0-3µm. La barre rouge correspond à la bande B11 infrarouge de Sentinel-2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.12 Exemple d'une image Formosat-2 du bassin de la Rheraya. La distinction neige/nuage peut se faire à l'oeil. . . . . . . . . . . . . . . . . . . . . . . . 18 2.11 Les différents produits de télédétection utilisés durant cette thèse. Formosat-2, Sentinel-2 et MODIS ont permis d'extraire les surfaces enneigées durant les saisons mentionnées. Tandis que Pléiades a permis de générer un modèle numérique terrain à très haute résolution spatiale. . . . . . . . . . . . . . . 34 2.12 Photographie de la station météorologique d'Oukaimeden. . . . . . . . . . 35 2.13 Évolution de la SCF en fonction du SWE avec les deux approches : a) Seuil fixe, b) SDC de Zaitchik and Rodell (2009) (Eq. 2.10). . . . . . . . . . . . 36 2.14 Principe de détermination de la matrice de confusion et du calcul de ses éléments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 Photograph taken near Oukaimeden in the Rheraya catchment on 14 December 2016. (31.194 o N, -7.584 o S). It illustrates the effect of slope and aspect on the snow cover area at decametric scales that is often observed in the High Atlas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

1.3. Suivi de l'évolution du manteau neigeux

Dans le contexte marocain compte-tenu des remarques précédentes, il semble judicieux d'opter pour une approche qui combine la télédétection et la modélisation distribuée afin

Ces méthodes se basent sur la minimisation d'une fonction de coût. Comme c'est difficile de la minimiser en utilisant le gradient quand les dimensions du vecteur de contrôle sont trop élevées, une technique numérique dite méthode est utilisée. Pour d'avantages informations, le lecteur peut se référer à[START_REF] Blayo | Advanced Data Assimilation for Geosciences : Lecture Notes of the Les Houches School of Physics : Special Issue[END_REF]

EnBS : Une méthode d'ensemble des filtres de Kalman (Ensemble Batch Smoother, en anglais)

c'est la première génération des données MERRA-2 que nous avons utilisées dans le chapitre 4

Equifinalité : la capacité des processus à atteindre le même état final, à partir de deux conditions initiales différentes

Remerciements

Acknowledgments This study was supported by the CNES TOSCA program. M. W. Baba has a doctoral grant funded by the CNES and the Région Occitanie (formerly Midi-Pyrénées). The Pléiades images acquisition was supported by the CNES through the ISIS program. The field work was funded by Agence Universitaire de la Francophonie through the NivAtlas project. We also acknowledge the Laboratoire Mixte TREMA for logistic and financial support. We thank Bastien Richard and Iñigo Irarrázaval for their participation to the GPS survey. We thank G.E. Liston for sharing the SnowModel code and his advices. The station data, the Formosat-2 snow maps and the GPS data used in this study are available under a Creative Commons Attribution Share-Alike 4.0. license. The Pléiades data are commercial data that can be obtained via http ://www.intelligence-airbusds.com/geostore/. Academic access exists through ISIS (http ://www.satelliteimageaccess.teledetection.fr/). Interested readers should contact the

Data

Despite of the importance of the snow cover and its high interannual variability, there is only one automatic snow depth record at the western edge of the catchment at Oukaimeden (3200 m.a.s.l) nearly on the water divide with the Rheraya catchment (Fig. 4.1). The snow depth was measured by using an acoustic sensor every half-hour and aggregated

corresponding author to request a specific agreement with CNES.We are grateful to Jessica Lundquist (editor), Nick Rutter (associate editor) and to the anonymous reviewers who provided many insightful comments that helped improve the quality of the paper.

Snowpack model

We used SnowModel [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF], a distributed energy-balance model, to simulate the snowpack evolution. The model was run at the hourly time step since the daily time step do not allow an accurate representation of the energy fluxes and determination of the precipitation phase [START_REF] Bair | Validating reconstruction of snow water equivalent in California's Sierra Nevada using measurements from the NASA Airborne Snow Observatory[END_REF]. SnowModel is composed of four sub-models : MicroMet, EnBal, SnowPack and SnowTran-3D. MicroMet distributes the meteorological forcing from station data on model grid using topography data [START_REF] Liston | A meteorological distribution system for highresolution terrestrial modeling (MicroMet)[END_REF]. EnBal solves the energy balance equation to compute the SWE evolution Liston and Elder (2006a). SnowPack computes the evolution of the snow depth and snowpack density Liston and Elder (2006a). SnowTran-3D computes the redistribution of the snow cover due to the wind [START_REF] Liston | Simulating complex snow distributions in windy environments using SnowTran-3D[END_REF].

In this work, SnowTran-3D was not activated because its high computational cost. However, the assimilation scheme is compatible with SnowTran-3D. In addition, the assimilation scheme is not dependent upon the snowpack model itself, i.e. submodels EnBal and SnowPack could be replaced by another implementation of a snowpack model. Although our assimilation scheme does not depend on the type of snowpack model, it is closely linked to MicroMet (as explained in Section 5.4.5). This is why we focus below on