# Machine-Learning-Based Downscaling of Hourly ERA5-Land Air Temperature over Mountainous Regions

Badr-Eddine Sebbar, Saïd Khabba, Olivier Merlin, Vincent Simonneaux, Chouaib El Hachimi, Mohamed Hakim Kharrou, Abdelghani Chehbouni

► To cite this version:

Badr-Eddine Sebbar, Saïd Khabba, Olivier Merlin, Vincent Simonneaux, Chouaib El Hachimi, et al. Machine-Learning-Based Downscaling of Hourly ERA5-Land Air Temperature over Mountainous Regions. Atmosphere, 2023, 14 (4), pp.610. 10.3390/atmos14040610. [https://doi.org/10.3390/atmos14040610](https://doi.org/10.3390/atmos14040610)

HAL Id: hal-04296516

https://hal.science/hal-04296516v1

Submitted on 20 Nov 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# Machine learning based downscaling of hourly ERA5-Land air temperature over mountainous regions

Badr-eddine Sebbar <fnref n="1" />,<fnref n="2" />* , Saïd Khabba <fnref n="1" />,<fnref n="3" />, Olivier Merlin <fnref n="2" />, Vincent Simonneaux <fnref n="2" />, Chouaib Elhachimi <fnref n="1" />, Mohamed Hakim Kharrou <fnref n="4" /> and Abdelghani Chehbouni <fnref n="1" />,<fnref n="2" />,<fnref n="4" />

**Abstract:** In mountainous regions, the scarcity of air temperature (Ta) measurements is a major limitation for hydrological and crop monitoring. An alternative to in-situ measurements could be to downscale the reanalysis Ta data provided at high-temporal resolution. However, the relatively coarse spatial resolution of these products (i.e., 9 km for ERA5-Land) is unlikely to be directly representative of actual local Ta patterns. To address this issue, this study presents a new spatial downscaling strategy of hourly ERA5-Land Ta data with a three-step procedure. First, the 9 km resolution ERA5 Ta is corrected at its original resolution by using a reference Ta derived from the elevation of the 9 km resolution grid and an in situ estimate over the area of the hourly Environmental Lapse Rate (ELR). Such a correction of 9 km resolution ERA5 Ta is trained using several machine learning techniques including Multiple Linear Regression (MLR), Support Vector Regression (SVR), and Extreme Gradient Boosting (Xgboost), and ancillary ERA5 data (daily mean, standard deviation, hourly ELR and grid elevation). Next, the trained correction algorithms are run to correct 9 km resolution ERA5 Ta, and the corrected ERA5 Ta data are used to derive an updated ELR over the area (without using in situ Ta measurements). Third, the updated hourly ELR is used to disaggregate 9 km resolution corrected ERA5 Ta data at the 30-meter resolution of SRTM's Digital Elevation Model (DEM). The effectiveness of this method is assessed across the northern part of the High Atlas Mountains in central Morocco through 1) k-fold cross-validation against five years (2016 to 2020) of in-situ hourly temperature readings, and 2) comparison to classical downscaling methods based on a constant ELR. Our results indicate a significant enhancement in the spatial distribution of hourly local Ta. By comparing our model, which included Xgboost SVR and MLR, to the constant ELR-based downscaling approach, we were able to decrease the regional root mean square error from approximately 3°C to 1.61°C, 1.75°C, and 1.8°C, reduce the mean bias error from -0.5°C to null, and increase the coefficient of determination from 0.88 to 0.97, 0.96 and 0.96 for Xgboost, SVR, and MLR respectively.

**Keywords:** Reanalysis; ERA5-Land; air temperature; downscaling; complex terrain; machine learning.

## 1. Introduction

Access to spatially and temporally consistent climate data at high spatial and temporal resolutions have progressively turned into a growing need in the 21st century, for being a paramount to numerous fields of study that investigate ecological, hydrological, and climate change processes, among others [1–7]. Using numerical weather models and data

**Citation:** Sebbar, B.; Khabba, S.; Merlin, O.; Simonneaux, V.; Elhachimi, C.; Kharrou M.H.; Chehbouni, A. Machine learning based downscaling of hourly ERA5-Land air temperature over mountainous regions. *Atmosphere* 2023, 1,0. https://doi.org/

Received:

Accepted:

Published:

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2023 by the authors. Submitted to *Atmosphere* for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Mohammed VI Polytechnic University, Center for Remote Sensing Applications, Ben Guerir, Morocco; badreddine.sebbar@um6p.ma; said.khabba@um6p.ma; chouaib.elhachimi@um6p.ma; abdelghani.chehbouni@um6p.ma</footnote>

<footnote n="2">² Centre d'Etudes Spatiales de la Biosphère (CESBIO), Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France; olivier.merlin@cesbio.cnes.fr; vincent.simonneaux@ird.fr</footnote>

<footnote n="3">³ LMFE, Physics department, Faculty of Sciences Semlalia, Cadi Ayyad University, Marrakech, Morocco; khabba@uca.ac.ma</footnote>

<footnote n="4">⁴ International Water Research Institute (IWRI), Mohammed VI Polytechnic University (UM6P), Ben Guerir, Morrocco; abdelghani.chehbouni@um6p.ma</footnote>

<footnote>* Correspondence: badreddine.sebbar@um6p.ma; Tel.: +212-653-01-12-60</footnote>
assimilation techniques to produce model-based reanalysis products is one viable strategy for generating climate datasets in light of this need [8–11]. Several international and local meteorological centers and data assimilation offices have collaborated over the past few decades to make numerous reanalysis products available to the public [12]. Examples of the most popular reanalysis products are: the ERA5 and ERA5-Land from European Centre for Medium Range Weather Forecasts (ECMWF); the second version of Modern-Era Retrospective Analysis for Research and Applications (MERRA2) [13] produced by NASA’s Global Modeling and Assimilation Office (GMAO); the second version of Climate Forecast System Reanalysis (CFSv2) from the National Centers for Environment Prediction and National Center for Atmospheric Research (NCEP/NCAR) [14]; NCEP/NCAR Global Reanalysis Products from NCEP and NCAR [15,16]; and the Japanese 55-year Reanalysis (JRA-55) [17] from the Japanese Meteorological Agency (JMA) [12,18]. The common key strength of the numerous existing reanalysis products resides in providing global data sets devoid of gaps, at high temporal resolution and over long time periods (generally over 3 or more decades). Still, reanalysis data frequently fail to simulate many of the processes that drive regional and local climate variability. Their limitations lie in their incapability of accurately depicting sub-km scale climate variables at the needed timescales and do not allow for proper representations of the local topography and sub-grid scale features that are essential in areas with complex terrain, microclimates or narrow mountain valleys, as highlighted by Holden et al. [19], Zhang et al. [20], Le Roux et al. [21], Alessi & DeGaetano [22], Zhang et al. [23]. When evaluated in contrast to observational data, the raw output data are regularly found to have systematic biases [24,25] limiting their usefulness for local applications [26]. There is consequently a need to make local-scale predictions more skillful by utilizing reanalysis data as input. In this context, a variety of techniques, such as downscaling methods, have been developed to bridge the gap between the scale at which data are available and the scale at which they are needed. The commonly used methods include dynamical downscaling and statistical downscaling [27].

Dynamical and statistical downscaling techniques are frequently used to refine coarser climate products to higher resolution [28,29]. The former is a widely used methodology to enhance the spatial information [30], in which a higher-resolution model, such as a regional climate model (RCM), can be driven by reanalysis data and run at spatial resolutions of up to a few meters’ projections (e.g. [31]), at which complex topography and smaller-scale processes are better represented [32]. This approach can give a very good simulation of local atmospheric conditions; however, it has significant computational cost [30,31,33,34]. Statistical downscaling methods, on the other hand, use statistical relationships to anticipate the evolution of local variables from large-scale variables. They are computationally less demanding and represent a more flexible alternative to dynamical downscaling. These methods have been shown to be effective in reproducing the fine-scale temperature variability over mountainous regions, particularly when using local observations (e.g., [1,35,36]).

This paper focuses on reanalysis air temperature (Ta) disaggregation over complex terrain since (i) it is one of the most important input variables in agro-environmental models and a crucial field for the vast majority of weather and climate applications, including climate change studies (e.g., [37,38]), and (ii) this variable is projected to change significantly in regions with irregular topography i.e., complex topography of mountain landscapes known to have a highly variable climate, with microclimates that can differ significantly from the surrounding area (e.g., [39,40]). Thus, having high-resolution Ta data over mountains allows for a better understanding of the complex microclimates that exist within mountain ranges, and can be particularly useful for predicting weather patterns and for understanding the impacts of climate change on these regions. Several studies describe the spatial interpolation methods used for downscaling in meteorology and climatology [37,41]. These techniques include nearest neighbor methods, splines, regression, kriging and cokriging but also machine learning techniques such as Artificial Neural Networks and Support Vector Machines [42–45]. None of these studies, however, focused on adjusting
reanalysis data to the regional real measured conditions prior to downscaling, nor worked on hourly time-step required for hydrological modeling relying on the availability of quality meteorological inputs at the simulation time step [46]. Recently, Sourp et al. [47] developed a snow reanalysis pipeline using downscaled ERA5 and ERA5-Land data. The downscaling is based on the MicroMet model [48,49] which performs spatial interpolation of meteorological variables using 100-m DEM [47]. Particularly, air temperature is downscaled to hourly timestep using the DEM and constant monthly Environmental Lapse Rates (ELRs).

Extending these previous ideas, a machine learning/statistical downscaling scheme is designed in this study to disaggregate hourly air temperature data with a 30-m spatial resolution from the 9-Km ERA5-Land Ta. The main originality relies on the assumption that the temporal variability of ELR should be taken into account for improving the spatial distribution of downscaled Ta estimates. The approach is tested in a steep-sided catchment in the western part of the High Atlas Mountains in central Morocco, where in situ Ta measurements are available from 2016 to 2020. The paper is organized as follows: the study area, data sets, and the methodology are presented in Section 2. Section 3 presents and discusses the results, while Section 4 outlines the principal conclusions.

## 2. Materials and Methods

### 2.1. Study area

The High Atlas is a large mountain range located in Morocco, stretching for 800 km in length and 60 km in width. It runs in a north-east to south-west direction and is known for its diverse range of elevations, from the lowest point of 1060-m above sea level to the highest peak in North Africa, Mount Toubkal, which reaches an elevation of 4167-m above sea level (Figure 1) [50,51]. The western part of the High Atlas is particularly notable for being a vital source of water for the northern plain of the Tensift catchment, specifically around the city of Marrakech [52]. The high-altitude regions of the mountain range are known for their low temperatures and sparse vegetation cover, with most agricultural activities concentrated along river valleys [53,54]. The Rheraya sub-basin (Figure 1), which is located 40 km south of Marrakech (between latitudes 30°05'N and 30°20'N, and longitudes 7°40'W and 8°00'W) and covers an area of 225 km², is one of the most intensely studied areas of the High Atlas Mountains. It represents a part of the Tensift Observatory in the frame of the SudMed [52] and the Joint International Laboratory LMI-TREMA [55] (https://www.lmi-trema.ma/) funded by the University Cadi Ayyad (UCA, Morocco) and the French Research Institute for Development (IRD, CESBIO Laboratory, France). It is considered to be globally representative of the western watershed of the High Atlas Mountains. The sub-basin contains four AWSs and is equipped with a variety of instruments to study the area.
<figure><img src="image_2.png" /><figcaption>Figure 1. Location of the Tensift basin and Rheraya sub-basin.</figcaption></figure>

## 2.2. Data set

### 2.2.1. Observed ground-based data

The air temperatures in the Rheraya sub-basin were measured on a semi-hourly basis from Automatic Weather Stations (AWSs) positioned throughout the sub-basin: Imskerbour (1404 m above sea level), Aremd (1940 m above sea level), Neltner (3207 m above sea level), and Oukaimden (3230 m above sea level). The temperature records for the period of 2016 to 2020 were converted from their original format to hourly time steps and any half-hour intervals with missing records from one or more AWSs were excluded. In order to ensure the accuracy of the data, the temperature records were checked for any excessive amounts of missing values, as outlined in the study of Dodson and Marks [56]. The missing values for the combined stations were ensured to not exceed 100 days per year. After the preprocessing step, the minimum number of hours kept per day for all years is 22 hours/day. The locations of the stations are illustrated in Figure 1 and Table 1 provides detailed information on the station names, heights, coordinates, yearly mean temperatures, number of observations, and frequency.

<table><caption>Table 1. Information regarding the four AWSs installed in the Rheraya sub-basin. The data collection period for all stations extends from January 1, 2016, to December 31, 2020.</caption><thead><tr><td>AWS</td><td>Latitude</td><td>Longitude</td><td>Elevation (m.a.s.l)</td><td>Tmean(°C)</td><td>No. of observations</td><td>Frequency</td></tr></thead><tbody><tr><td>Imskerbour</td><td>31.21018°</td><td>-7.93972°</td><td>1404</td><td>15.06</td><td>40870</td><td>30 min</td></tr><tr><td>Aremd</td><td>31.12948c</td><td>-7.91967°</td><td>1940</td><td>12.1</td><td>43848</td><td>30 min</td></tr><tr><td>Neltner</td><td>31.06579°</td><td>-7.91389°</td><td>3207</td><td>6.04</td><td>43829</td><td>30 min</td></tr><tr><td>Oukaimden</td><td>31.19328°</td><td>-7.86546°</td><td>3230</td><td>5.85</td><td>42644</td><td>30 min</td></tr></tbody></table>

### 2.2.2. Reanalysis data

For this study, the most advanced global reanalysis data produced in Europe, specifically optimized for land surface applications, was used. The dataset used is the ERA5-Land enhanced global dataset for the land component of the fifth generation of European Reanalysis, which is freely available on the website https://cds.climate.copernicus.eu [57]. In
comparison to ERA5 and older ERA-Interim, the ERA5-Land dataset has the advantage of enhanced horizontal resolution of 9-km (released on a regular 0.1° x 0.1° grid) compared to 31-km (ERA5) and 80-km (ERA-Interim), while maintaining the same hourly temporal resolution as ERA5 [18]. The high temporal and spatial resolutions of ERA5-Land, and the consistency of the fields produced, make it a valuable dataset for diverse applications related to water resources, land and environmental management. The variable of interest in this study is the air temperature at 2-m above the surface of land, sea or in-land waters, which is calculated through interpolation between the lowest model level and the Earth's surface, considering atmospheric conditions. Hourly ERA5-Land temperature data was downloaded and processed to be consistent with the measured data screening. Additionally, ERA5-Land temperature's daily mean, minimum, maximum and standard deviation were computed as ancillary data for the entire study period.

### 2.2.3. Digital Elevation model

To achieve fine-scale disaggregation, we used the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global digital elevation model (DEM) with 30-m resolution (https://earthexplorer.usgs.gov/). The related SRTM 1 Arc-Second tile (SRTM1N31W008V3) was used to generate a DEM subset of the Rheraya sub-basin (figure 1).

## 2.3. Methodology

In this section, we outline the process for enhancing the spatiotemporal downscaling of Ta. We first explain the use of machine learning models to correct ERA5-Land Ta (hereafter referred to as Ta_5) using in-situ hourly Ta and ELR readings (Ta_st and ELR_st), resulting in corrected Ta_5 (Ta_5_corr). Then, we describe the process of using Ta_5_corr to downscale temperatures at a 30-meter resolution using a DEM, producing Ta_disagg_ML. This final product will be validated against five years of in-situ hourly Ta_st readings and compared to two other downscaling methods (Annual ELR average and MicroMet model) to evaluate the improvement made. Further details on each step are provided in subsequent subsections.

### 2.3.1. 1st step: Ta_5 correction

The process of correcting Ta_5 starts with the creation of a reference Ta (Ta_5_ref) corresponding to each 9-km ERA5-Land grid elevation, utilizing only ground data, specifically, hourly measured Ta_st and ELR_st. The Ta_5_ref is aligned with the measured Ta_st and ELR and is intended to be more accurate than the one provided by ERA5-Land, serving as the target to be achieved prior to downscaling. This step is illustrated in figure 2. In the second step, using only ERA5-Land Ta_5, a set of variables is derived, and which may be correlated with the local disaggregated temperature that is intended to be produced. This set of variables will then be utilized in the machine learning approaches. In the third step, Ta_5 is corrected to match the Ta_5_ref (9-km spatial resolution) using machine learning models. In these models, the estimated value is Ta_5_corr, Ta_5_ref is the dependent variable, and the independent variables include Ta_5 and the selected variables from step 2.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example of Ta_5_ref estimates for ERA5-Land grid elevations based on observed hourly ELR_st (slope). The dashed black line represents the regression line of measured temperature to elevation. The red dashed lines show the difference of ERA5-Land Ta_5 to Ta_5_ref (what it should be).</figcaption></figure>

The Ta_st measured by AWSs are plotted against their corresponding elevations, and linear regressions are used to calculate the slope hourly ELR_st and the intercept b_st (representing air temperature at sea-level). These values are then used to interpolate hourly Ta_5_ref for elevations of ERA5-Land grid points (9-km spatial resolution) over the period of interest (from 2016 to 2020). The equation that governs this interpolation is as follows (equation 1):

$$Ta\_5\_ref = ELR\_st \times E5 + b\_st \quad (1)$$

E5 being the elevation of ERA5-Land grid point in meters.

These Ta_5_ref values are then used to calibrate machine learning models as a dependent variable to correct Ta_5. The input variables include Ta_5 and a set of variables potentially correlated to the local disaggregated temperature that is intended to be produced. The input features for predicting a specific variable may be highly correlated with one another, resulting in a processing and computational time loss. In addition, those input features may not always be correlated with the target variable, which can result in an overfitting of the constructed model. In other words, the learned model would be a better fit for training data than test data [58]. To avoid those problems, carrying out a correlation analysis holds the key to decide inputs to keep or to exclude. The Pearson Correlation Coefficient (PCC) can be used to calculate the correlation between candidate input variables Xi and targeted variable Y. The PCC is by definition the covariance of X_i and Y over the product of their standard deviations σ_X_i and σ_Y. It ranges from -1 to +1, where a value of -1/+1 implies that X_i is completely negatively/positively linearly correlated to Y, and a value of 0 indicates absolute absence of correlation between the two variables. In most cases, a high absolute value of PCC (often greater than 0.8) indicates strong correlation [58]. The expression of PCC is given in equation 2:

$$PCC = \frac{COV(X_i, Y)}{\sigma_{X_i} \times \sigma_{Y}} \quad (2)$$

The candidate variables selected for conducting the correlation analysis are hourly Ta_5, hourly ELR calculated using Ta_5 (hereafter referred as ELR_E5), daily Ta_5 means, minimums, maximums, standard deviations, and ERA5-Land grid points elevations. These variables will be used to examine the correlation with the targeted variable, Ta_5_ref. The goal is to use the selected input variables, all of which are sourced from ERA5-Land data, to predict a more accurate corrected Ta_5_corr which will be applied for downscaling. We have chosen to test three different models for predicting Ta_5_ref: (1) a basic multiple input
linear regression method known as MLR, (2) the popular and widely used SVR model, and (3) one of the newest machine learning methods the Xgboost algorithm, which is known for its exceptional predictive abilities. Next is a brief theoretical explanation of the operation and functioning of the models.

* MLR

In MLR, multiple independent variables are used to describe the behavior of the dependent variable [59]. It is an extension of simple linear regression, and it describes the relationship between two or more explanatory variables and a response variable by fitting a linear equation to observed data. Each value of the independent variable corresponds to a prediction value for the dependent variable. A good MLR model should be able to explain a majority of the variance in the dependent variable with the smallest number of independent variables possible. For a more detailed explanation of MLR theory, the reader is encouraged to refer to the work of Helsel and Hirsch [60].

* SVR

SVR is a branch of Support Vector Machine (SVM) that is widely used as a regression technique (detailed description of SVM can be found in several works (e.g., [61–63]). SVR finds a multivariate regression function that predicts a desired output property or dependent variable Y based on a set of input independent variables X (NxM) and Y (M). The main difference between SVR and MLR is that in SVR, the original input space (which is usually non-linearly related to the targeted variable) is mapped onto a higher dimensional feature space using a kernel function (such as Linear, Radial Basis Function, Polynomial, sigmoid) to find an optimal hyperplane to separate the sample points. The full description of SVR equations is not included here but can be found in works such as [64], [65], and [66].

* Xgboost

Xgboost, proposed by Chen et al. in 2015 [67], is an alternative method for predicting a response variable based on certain covariates. It is similar to the well-known Random Forest method, it builds classification and regression trees one by one, but instead of making a decision based on a final vote, each subsequent model (tree or base learner) is trained using the mistakes of the previous one. This technique is becoming increasingly popular due to its design and ability to speed up training time using various techniques such as parallel computing and sparsity-aware split-finding. For more details the reader is referred to the following [34,67,68]

All the previously mentioned algorithms were implemented using the Python library "Scikit-learn" developed by Pedregosa et al., in 2011 [69]. Scaling was performed prior to using SVR kernel methods as they are based on distance, this was done to facilitate learning and prevent features with the largest range from dominating the computations. The "RobustScaler" method was used as it can handle outliers. The performance of the machine learning models heavily depends on the hyperparameter values, therefore, a significant step was determining the optimal values for the model through hyperparameter tuning. This was done using the "Scikit-learn" library's Grid Search function, which considers multiple hyperparameter combinations and chooses the one that returns the lowest error score. Since MLR model does not have any hyperparameters to tune, only SVR's and Xgboost's hyperparameters were tuned. The Grid search function also includes a pre-defined k-fold cross validation method [70–74], where each fold serves as a single hold-out test fold, and the model is built using the remaining k-1 folds. Grid search methodology with 5-fold cross-validation was applied to get the optimal model parameters for SVR and Xgboost, meaning that during the cross-validation process, 4 years of data were used for calibration and 1 year of data for validation.

### 2.3.2. 2nd step: Disaggregation

Climate impact studies frequently use a constant Ta lapse rate at specific locations, that we hereafter refer to as ELR_cst and which is equal to -6.5 °C.km⁻¹ (e.g., [75–78]). However, this rate can vary significantly depending on factors such as location, season, and
time of day. Studies have shown that the temperature ELR can range from -9.8 °C.km⁻¹ to -10 °C.km⁻¹ in dry conditions (the dry adiabatic lapse rate), and values that are shallower or equal to -6.5 °C.km⁻¹ generally represent moist adiabatic conditions [28,79–81]. This variability have been measured in our study area as it is shown in figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3. Pronounced ELR's hourly temporal variability, measured using AWSs record over the period of interest (from January 1st, 2016, to December 31st, 2020).</figcaption></figure>

To account for this variability, the correction of ERA5-Land Ta_5 on an hourly basis enables tracking actual local ELR values. Once the models predict corrected temperature values Ta_5_corr, new hourly temperature lapse rates (ELR_corr) are computed through linear regression, and then the Ta_5_corr (9-km) is downscaled to the DEM of the area of interest (30-m) using those corrected values instead of the original ones. The equation used for this downscaling process is shown in Equation 4 and the classic constant Ta lapse rate method's formula displayed in equation 3.

$$Ta\_disagg\_cst = Ta\_5 + ELR\_cst \times (DEM - E5) \quad (3)$$

$$Ta\_disagg\_ML = Ta\_5\_corr + ELR\_corr \times (DEM - E5) \quad (4)$$

In addition to using machine learning and constant Ta lapse rate approach to down-scale Ta_5, the MicroMet model was also applied for comparison. The MicroMet model is a high-resolution meteorological distribution model designed to produce high-resolution meteorological data such as air temperature, humidity, wind, radiation, and precipitation for use in running spatially distributed terrestrial models over a variety of landscapes. It uses established relationships between meteorological variables and the surrounding landscape to distribute those variables in a computationally efficient and physically plausible way. Specifically for air temperature, the MicroMet model first adjusts the Ta_5 values to sea level using the formula (equation 5):

$$Ta\_0 = Ta\_5 + ELR\_month \times (E5 - E0) \quad (5)$$

Ta_0 and ELR_month being the Ta adjusted to sea-level and the monthly values of the ELR respectively (see Table 2), where the ELR_month values vary depending on the month of the year [82] or are calculated based on data from nearby stations. The sea-level Ta_0 values are then interpolated to the model grid using the Barnes objective analysis method [83]. The gridded topography data and ELR_month are then utilized to adjust the sea-level gridded temperatures to the elevations provided by the DEM, using the equation provided in equation 6:

$$Ta\_disagg\_MM = Ta\_0 + ELR\_month \times (DEM - E0) \quad (6)$$
<table><caption>**Table 2.** Air temperature ELR (°C.km⁻¹) variations, for each month of the year, in the Northern Hemisphere [82].</caption><thead><tr><th>Month</th><th>Jan.</th><th>Feb.</th><th>Mar.</th><th>Apr.</th><th>May</th><th>Jun.</th><th>Jul.</th><th>Aug.</th><th>Sep.</th><th>Oct.</th><th>Nov.</th><th>Dec.</th></tr></thead><tbody><tr><td>ELR</td><td>4.4</td><td>5.9</td><td>7.1</td><td>7.8</td><td>8.1</td><td>8.2</td><td>8.1</td><td>8.1</td><td>7.7</td><td>6.8</td><td>5.5</td><td>4.7</td></tr></tbody></table>

### 2.3.3. 3rd step: Validation and results assessment

The quality of the final products (i.e., the downscaled Ta_disagg_ML) was evaluated through ins-situ validation and comparison to the other two scenarios (Ta_disagg_MM and Ta_disagg_cst) using statistical parameters. Three simulation evaluation scores were used: Root Mean Square (RMSE), coefficient of determination (R²) which is the square of the previously described PCC (Pearson's Correlation Coefficient), and the Mean Bias Error (MBE) [84]. The scores were computed for each AWS for validation. The mathematical expressions of the above scores are presented in equations 7 and 8 (R² is the square of PCC in equation 2).

$$ MBE = \frac{1}{N} \sum_{i=1}^{N} (Pred_i - Obs_i) \qquad (7) $$

$$ RMSE = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (Pred_i - Obs_i)^2} \qquad (8) $$

Predᵢ being the predicted value and Obsᵢ the measured one. The above-mentioned steps and methodology description have been summarized in the following flowchart (figure 4.). It provides a clear and concise summary of the method and can serve as a guide to understand and replicate the study's methodology.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Flowchart of the methodological approach.</figcaption></figure>
# 3. Results

302

## 3.1. Ta_5 correction

303

* Reference temperature Ta_5_ref

304

The obtained hourly reference Ta_5_ref values were compared to the Ta_5 values sourced directly from ERA5-Land data. The comparison was done over all ERA5-Land grids and the entire study period (2016 to 2020). The results of this comparison are shown in figure 5. The mean R², RMSE, and MBE of Ta_5_ref and Ta_5 were found to be 0.88, 2.51 °C, and -0.48 °C, respectively. The results indicate that the predicted values closely follow the reference values, however, the difference between the two can reach up to approximately 10 °C.

305

306

307

308

309

310

311

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of ERA5-Land's original Ta_5 and reference Ta_5_ref air temperatures.</figcaption></figure>

The next figure ( figure 6) illustrate a temporal comparison of Ta_5_ref to the original Ta_5 for two ERA5-Land grids. The lines are plotted on top of each other and the difference between the two temperatures can be easily observed. The figure shows an example of the comparison for two ERA5-Land grids over the first two and a half months of 2016, and similar behavior is observed throughout the study period.

312

313

314

315

316
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Hourly Comparison of Ta_5 and reference Ta_5_ref over time (as dashed red line and black line respectively). Example of the ERA5-Land grid situated at: (a) 7.9°W and 31.1°N, and (b) 7.9°W and 31.2°N.</figcaption></figure>

The plot indicates that the trend of the two variables is similar, meaning that they both increase or decrease at the same rate over time. However, the amplitude of the Ta_5_ref variable is less than the amplitude of the original reanalysis Ta_5, meaning that the range of temperatures it covers is smaller. This suggests that the reanalysis Ta_5 has higher amplitude of Ta variations than what it should be over the study area. The corrections to be applied to the Ta_5 are then to adjust for the bias that may present in the reanalysis dataset. This bias can be caused by errors in the input data, topographical effects, the modeling approach or in the assimilation of observations. The bias can also be caused by the lack of representation of the complex topography or urbanization of the study area in the reanalysis dataset. At first, we attempted to debias/correct the Ta_5 using simple linear regression, modeling daily temperature changes as a sinusoidal function, and constant (positive or negative) bias correction prior to downscaling. However, these methods did not yield significant improvement and were not practical for the study area, hence the choice of the machine learning approaches. As we stated in the methodology section a correlation analysis was carried out to select proper input variables prior to Ta_5 correction, and thus based only on a set of variables independent from in-situ data (Ta_5 and its derivates, as well as ERA5-land grid points elevations).

* Correlation analysis and feature selection

Figure 7 depicts the results of the correlation analysis. To test for correlated input variables, the independent targeted variable Ta_5_ref was also introduced to the correlation matrix. The latter shows that Ta_5, Ta_5's daily minimum, Ta_5's daily maximum, and Ta_5's daily mean are all highly correlated, with PCCs of 0.95, 0.87, 0.88, and 0.89, respectively. Moderately low to low correlations are found for standard deviation (PCC= 0.44), ELR_E5 (PCC= -0.18), and E5 (ERA5-land's grid elevation) (PCC= -0.26). Nonetheless, given our emphasis on finer resolution and higher precision, keeping those inputs appears to be very appropriate as long as they are not very close to null (under 0.05 for instance).
<figure><img src="image_8.png" /></figure>

**Figure 7.** Correlation matrix results. PCC value of each two variables is shown in the boxes corresponding to their "coordinates". ELR_E5, Std and E5 being the hourly ELR issued from Ta_5, the daily standard deviation, and ERA5-land grids elevation respectively.

The correlation matrix also shows that the daily mean of Ta_5 have almost perfect correlation with both the daily minimum and maximum values, with correlation coefficients of 0.97 and 0.99 respectively. Additionally, among the three, the daily mean showed the best correlation to the targeted variable Ta_5_ref (correlation coefficient of 0.89), thus only the mean was kept. The final set of retained input variables for predicting Ta_5_ref values (i.e., correcting Ta_5) consisted of five variables: Ta_5, Ta_5's daily mean and standard deviation, hourly ELR_E5 and E5 (ERA5-land's grid elevation). It is worth noting that while elevation remains constant over time, it varies from one ERA5-Land grid to the next, hence its inclusion was entirely justified.

* Machine learning outcome

The three scatterplots of figure 8 compare the predictions of temperature made by the tested machine learning algorithms, MLR, SVR, and Xgboost, with the reference Ta_5_ref. Overall, the results show a good level of agreement between the predictions of the three models and the targeted reference Ta_5_ref.
<figure><img src="image_9.png" /><figcaption>Figure 8. Comparison of Machine Learning predictions for Ta_5_ref temperature</figcaption></figure>

The MLR-based model had a RMSE of 1.34 °C, a R² of 0.97, and a quasi-null MBE. The fitting parameters of the MLR model are the coefficients of the regression equation used to predict the reference Ta_5_ref. They represent the contribution of each input feature in the linear equation. The specific values found for these fitting parameters are as follows: 0.507 for hourly Ta_5, 0.477 for daily mean, -3.17×10⁻³ for ERA5-Land grids elevation, -186.71 for hourly ELR_E5, -0.329 for daily standard deviation, and 6.27 for the intercept.

The SVR model used the Radial Basis Function (RBF) kernel which is known to provide good general performance as reported in previous studies such as Zaidi (2015) and Parveen et al. (2016). The grid search methodology along with 5-fold cross validation was utilized to find the best values for the SVR model parameters, such as C, ε, and γ. A wide range of permutations were tried and tested, such as C [2⁻², 2¹²], γ [2⁻¹², 2²] and ε [2⁻¹², 2⁴]. The statistical evaluation mean parameters for the best fitted SVR model were found to be RMSE = 1.20 °C, R² = 0.97, and MBE=0 °C using the Python package scikit-learn and the rules of "Lesser is better" for the RMSE and MBE and "Greater is better" for R². The best parameters found were C = 1, γ = 'scale' and ε = 0.02.

The grid search methodology was also applied to the Xgboost algorithm to find the best evaluation metrics (lowest RMSE and MBE and highest R²). An analysis of Aarshay's (2016) work was used as a reference to determine typical values of learning rate, maximum depth, minimum child weight, gamma, subsample, and colsample by tree, such as: [0.01,0.2], [3,10], [1,6], [0.1,0.2], [0.5,0.9], and [0.5,0.9]. The best fit was found when the following settings were used: learning rate= 0.4, maximum depth= 6, minimum child weight= 1, subsample= 1, colsample by tree= 1, and a "number of estimators" of 2000. The results from the Xgboost model were superior to those from the SVR model and MLR model, with an RMSE of 0.83 °C, R² of 0.99 and MBE of 0 °C respectively. Table 3 displays the specific outcomes for the three scoring parameters from the various cross-validation folds. Overall, we see a consistent pattern of model behavior throughout each fold change process, indicating that the models are well calibrated and are not overfitting.
<table><caption>**Table 3.** Detailed corss-validation results.</caption><thead><tr><th rowspan="2">Cross-validation (years)</th><th colspan="3">MLR</th><th colspan="3">SVR</th><th colspan="3">Xgboost</th></tr><tr><th>RMSE</th><th>R<sup>2</sup></th><th>MBE</th><th>RMSE</th><th>R<sup>2</sup></th><th>MBE</th><th>RMSE</th><th>R<sup>2</sup></th><th>MBE</th></tr></thead><tbody><tr><td>2016</td><td>1.3878</td><td>0.9500</td><td>0.3740</td><td>1.2213</td><td>0.9690</td><td>0.2622</td><td>0.8411</td><td>0.9877</td><td>0.0240</td></tr><tr><td>2017</td><td>1.3542</td><td>0.9584</td><td>-0.1944</td><td>1.2456</td><td>0.9654</td><td>-0.2187</td><td>0.8232</td><td>0.9874</td><td>-0.0123</td></tr><tr><td>2018</td><td>1.2935</td><td>0.9670</td><td>-0.0129</td><td>1.1828</td><td>0.9733</td><td>-0.0112</td><td>0.7891</td><td>0.9870</td><td>-0.0116</td></tr><tr><td>2019</td><td>1.3310</td><td>0.9654</td><td>0.1079</td><td>1.2174</td><td>0.9696</td><td>0.1109</td><td>0.8260</td><td>0.9872</td><td>-0.0031</td></tr><tr><td>2020</td><td>1.3234</td><td>0.9660</td><td>-0.2734</td><td>1.1789</td><td>0.9750</td><td>-0.1608</td><td>0.8139</td><td>0.9871</td><td>-0.0104</td></tr></tbody><tfoot><tr><td>Mean</td><td><strong>1.34</strong></td><td><strong>0.97</strong></td><td><strong>0.002</strong></td><td><strong>1.21</strong></td><td><strong>0.97</strong></td><td><strong>-0.004</strong></td><td><strong>0.83</strong></td><td><strong>0.99</strong></td><td><strong>-0.003</strong></td></tr></tfoot></table>

To sum up, exceptional Ta_5 correction performance of the Xgboost model in predicting the reference Ta_5_ref was observed. The high R² value and low RMSE and MBE values indicate a better fit compared to the MLR and SVR models. Additionally, the Xgboost model stands out for its combination of both speed and accuracy which is a significant advantage.

## 3.2. Ta_5_corr downscaling

In this section, we present the results of our study on downscaling the 9-km ERA5-land's Ta_5 using three different scenarios. As a reminder, the three scenarios explored are our own method, the machine learning-based method, and a comparison to classic downscaling approaches, the MicroMet model, as well as the often-used constant ELR method (ELR_cst). As previously stated, the machine learning method was used to correct the Ta_5 values, and new values for ELR were calculated from the corrected Ta_5_corr values. These ELR_corr values are then used for the downscaling of the latter temperature. The results of each scenario will be discussed in detail, and the comparison between them will be highlighted.

Figure 9 highlight improvements made on ELR_corr estimations posterior to Ta_5 correction. The first subplot is a scatterplot of the ELR issued from non-corrected Ta_5 against the measured ELR_st from AWSs. The second subplot is a scatterplot of the machine learning-based corrected ELR_corr against the measured ELR_st from AWSs (we are only showcasing the ELR-corr based on Ta_5_corr corrections made using Xgboost model, as it had the most favorable outcome).

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of non-corrected and corrected Ta_5 resulting ELRs.</figcaption></figure>

The scatterplots show that there is a significant improvement in the agreement with the measured ELR_st when using the machine learning-based approach ELR_corr. The R² value for the ELR_corr is 0.78, which is significantly higher than the R² value of 0.41 for the non-corrected ELR_E5. This indicates that the ELR_corr model has a better ability to accurately predict the measured ELR_st values from the AWSs. Constant ELR_cst and Micromet model's monthly values ELR_month were not compared as we would only
get horizontal lines given the fact that they are constant and the measured value ELR_st exhibits huge spatial temporal variability.

The results presented in Figure 10 show the performance of the three followed approaches for downscaling temperatures at 30-m resolution. The scatterplots compare the downscaled temperatures from each approach to the measured validation temperatures (Ta_st) from the four AWSs Imeskerbour, Aremd, Neltner, and Oukaimden. The AWSs are displayed in columns, while the rows indicate the approach followed.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Performance evaluation of machine learning-based ERA5-Land's temperature downscaling against traditional Methods using In-Situ measurements.</figcaption></figure>

The first approach, using ELR_cst the constant elevation-based lapse rate and Ta_5 the original ERA5-land's temperature data, performed poorly, as expected, yielding a RMSE of 3.11 °C, a coefficient of determination R² of 0.81, and a MBE of -0.55 °C. Using the MicroMet model, the second approach did not improve the predictions either although it outperforms the constant lapse rate model, with overall performance estimates of 2.71 °C, 0.85, and -0.40 °C for RMSE, R², and MBE, respectively.

The new machine learning based approach, which corrects Ta_5 temperature and lapse rate data prior to downscaling (Ta_5_corr and ELR_corr), showed a satisfying improvement in the match between downscaled and measured temperatures. The intercomparison of the three machine learning models (Xgboost, SVR, and MLR) revealed that the Xgboost model had the best performance, with a RMSE of 1.61 °C, a R² of 0.97, and a MBE of 0 °C. The SVR model had a slightly worse performance with an RMSE of 1.75 °C, a R² of 0.96 and a MBE of 0 °C, but it took significantly more time to compute. The MLR model had the lowest performance with a RMSE=1.8 °C, R² =0.96, and MBE=0 °C, but still presents a
satisfying improvement compared to constant elevation-based lapse rate and MicroMet models. The next table (table 4), provides further details on the downscaling performance metrics by station and approach. The table shows that overall, the constant lapse rate elevation-based approach and the MicroMet model present consistent RMSE for all the stations, however, the MBEs differ. These differences in relation to the measurements can be considered quite important, especially if the downscaled product is intended to be used as input for fine-scale models.

<table><caption>**Table 4.** Downscaling performance metrics by station and approach.</caption><thead><tr><th rowspan="3">AWS</th><th rowspan="3">Cst<br>ELR</th><th colspan="4">RMSE</th><th colspan="4">R<sup>2</sup></th><th colspan="4">MBE</th></tr><tr><th colspan="3">MicroMet</th><th rowspan="2">Cst<br>ELR</th><th colspan="3">MicroMet</th><th rowspan="2">Cst<br>ELR</th><th colspan="3">MicroMet</th><th rowspan="2">Cst<br>ELR</th></tr><tr><th>MLR</th><th>SVR</th><th>Xgboost</th><th>MLR</th><th>SVR</th><th>Xgboost</th><th>MLR</th><th>SVR</th><th>Xgboost</th></tr></thead><tbody><tr><td>Imskerbour</td><td>2.45</td><td>2.71</td><td>1.82</td><td>1.86</td><td>1.75</td><td>0.90</td><td>0.88</td><td>0.95</td><td>0.94</td><td>0.95</td><td>-0.43</td><td>0.28</td><td>0.34</td><td>0.42</td><td>0.34</td></tr><tr><td>Aremd</td><td>3.09</td><td>2.47</td><td>2.05</td><td>1.95</td><td>1.77</td><td>0.84</td><td>0.9</td><td>0.93</td><td>0.94</td><td>0.95</td><td>-2.14</td><td>-0.70</td><td>-0.50</td><td>-0.45</td><td>-0.49</td></tr><tr><td>Neltner</td><td>3.29</td><td>3.00</td><td>1.61</td><td>1.55</td><td>1.41</td><td>0.76</td><td>0.81</td><td>0.94</td><td>0.95</td><td>0.95</td><td>-0.49</td><td>-0.67</td><td>0.10</td><td>0.02</td><td>0.08</td></tr><tr><td>Oukaimden</td><td>3.47</td><td>2.67</td><td>1.68</td><td>1.62</td><td>1.47</td><td>0.75</td><td>0.82</td><td>0.94</td><td>0.94</td><td>0.95</td><td>0.86</td><td>-0.54</td><td>0.10</td><td>0.00</td><td>0.08</td></tr></tbody></table>

On the other hand, it is noted that all metrics are improved for all stations when using machine learning approaches. Additionally, it can be observed that the metrics for higher elevations (Oukaimden and Neltner) are better than those in lower elevations (Imskerbour and Aremd). This could be explained by several factors such as the larger differences in temperature between the high and low elevations, or a better alignment to regression lines, and hence better corrected Ta_5_corr values. It could also be due to the fact that the machine learning models are able to capture the complex interactions between temperature and ERA5-land grids elevation in these regions more effectively.

In conclusion, the results of this study indicate that the present machine learning-based downscaling technique has great potential for disaggregating ERA5-Land Ta_5 coarse 9-km resolution to the DEM's 30-meter resolution, particularly in harsh and difficult-to-access mountainous regions. The use of machine learning models improved the performance of the downscaling process and the match between predicted and measured Ta. This approach outperforms the traditional constant elevation-based lapse rate and MicroMet model. Additionally, the Xgboost model was found to be the best option for reproducing this methodological approach as it performed better and faster than the other two models (MLR and SVR).

The illustration presented in the next figure (figure 11) depicts an example of mapping across the study region and summarizes the strategy followed to create a high-resolution downscaled air temperature based on 9-km ERA5-Land's Ta_5 maps once the models are calibrated.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** High resolution temperature mapping of mountainous regions using machine learning-based downscaling of ERA5-Land's T2m data. Example showing Xgboost in action across the Rheraya basin on October 10, 2021, at 11 a.m., after the initial calibration of the model.</figcaption></figure>
## 4. Discussion

The correction of ERA5-Land Ta_5 data through the application of machine learning techniques resulted in an enhanced spatial distribution of downscaled Ta estimates. The improvement was demonstrated through comparison with two classical downscaling methods: the annual average and the MicroMet model. To summarize, the process began with the creation of Ta_5_ref temperatures calculated for each ERA5-Land grid points elevation to match the observed local temperature-elevation relationship. In simpler terms, Ta_5_ref is a 9-km adjusted version of the ERA5-Land Ta_5 and a more accurate representation of the actual measurements of Ta_st and ELR_st. Hence, Ta_5_ref served as the desired outcome for the correction of Ta_5.

The gap between Ta_5 and Ta_5_ref values was filled through machine learning. Three different machine learning techniques, MLR (simple), SVR (relatively complex), and Xgboost (recent), were selected to make the prediction of Ta_5_ref. A correlation analysis was performed to determine the input variables that could be correlated to Ta_5_ref. These candidate input variables were all derived from the ERA5-Land data, meaning that once the models were calibrated, the Ta_5 temperature was corrected using its own data to align with the observed local temperature-elevation relationship before downscaling. The results of the correlation analysis showed that the set of input variables is included in addition to hourly Ta_5: hourly ELR_E5, the mean and standard deviation of daily Ta_5, and the elevation of the ERA5-Land grid points. The predicted/corrected values at 9-km spatial resolution, referred to as Ta_5_corr, were validated against Ta_5_ref and showed significant improvement. The original gap between Ta_5 and Ta_5_ref was quantified as having a RMSE of 2.51°C, R² of 0.88, and MBE of -0.48°C. The MLR-based model showed a correction with a RMSE of 1.34°C, R² of 0.97, and a near-zero MBE. The best fit SVR model had a RMSE of 1.20°C, R² of 0.97, and MBE of 0°C. The Xgboost model performed even better, with an RMSE of 0.83°C, R² of 0.99, and MBE of 0°C, surpassing the results from the SVR and MLR models. The Ta_5_corr values at 9-km spatial resolution, more aligned with local measurements than the original Ta_5, were then used to calculate ELR_corr values. The resulting ELR_corr values were plotted against measurements and showed a R² of 0.78 and a RMSE of 0.001°C/km. The final product, the disaggregated Ta_5_disagg, was obtained by using Ta_5_corr and ELR_corr in conjunction with a 30-m DEM.

The downscaling results showed a satisfying improvement in the match between downscaled Ta_5_disagg and measured Ta_st. The intercomparison of the three machine learning models (Xgboost, SVR, and MLR) revealed that the Xgboost model had the best performance, with a RMSE of 1.61 °C, a R² of 0.97, and a MBE of 0 °C. The SVR model had a slightly worse performance with an RMSE of 1.75 °C, a R² of 0.96 and a MBE of 0 °C, but it took significantly more time to compute. The MLR model had the lowest performance with a RMSE=1.8 °C, R² =0.96, and MBE=0 °C, but still presents a satisfying improvement compared to constant elevation-based lapse rate and MicroMet models. These differences in relation to the measurements can be considered quite important, especially if the downscaled product is intended to be used as input for fine-scale models.

The limitation of this method is that it needs a starting point, i.e., the machine learning models must be first calibrated accordingly to the reference temperature Ta_5_ref that is calculated through in-situ measurements and ERA5-Land grid points elevations. The primary benefit, however, is that this is one of the few works that successfully downscales ERA5-Land Ta_5 to an hourly time-step, is applicable throughout all seasons, and captures both diurnal and regional temperature fluctuations. Moreover, once the models are calibrated over a specific area, they can be used independently of any knowledge of in-situ measurements as wa previously mentioned, the inputs consist solely of ERA5-Land Ta_5, its derived products (hourly ELR, daily mean and standard deviation), in addition to ERA5-land grids points elevations.
## 5. Conclusions

The ERA5-Land Ta_5 data has been improved through the use of machine learning techniques in downscaling. The correction process started with the creation of Ta_5_ref, which is a 9-km adjusted version of the ERA5-Land Ta_5, better representing the actual temperature measurements. The gap between Ta_5 and Ta_5_ref was filled through machine learning using three models: MLR, SVR, and Xgboost. The results showed that the Xgboost model performed the best, surpassing the SVR and MLR models. The downscaled product showed significant improvement compared to the one obtained through classic downscaling approaches (constant ELR and MicroMet model). The primary benefit of this method is that it can accurately downscale to an hourly time-step, is applicable throughout all seasons, and captures diurnal and regional temperature fluctuations. However, the models must be calibrated for a specific area before use. Overall, this method presents a promising solution for improving the accuracy of temperature data downscaling and can be used for other climate studies.

In perspective, assessment of the added value of this novel machine learning based method for hydrological applications is considered (e.g., reference evapotranspiration over mountains). Another avenue would be the extension of the use of machine learning models to downscale other meteorological variables (e.g., wind speed, relative humidity...). Finally, although the time-window would be more restrained, we can also consider the use of machine learning based methods on ERA5-Land's Land Surface Temperature (LST) to reproduce high resolution satellite product such as thermal based Landsat-8 LST.

**Author Contributions:** All the authors (S.B. and K.S. and M.O. and V.S. and EH.C. and C.A) have contributed substantially to this manuscript. "Conceptualization, S.B.; methodology, S.B.; software, S.B. and EH.C.; validation, S.B.; formal analysis, S.B. and K.S. and M.O. and V.S. and M.H.K and C.A.; investigation, S.B. and K.S. and M.O. and V.S. and M.H.K and C.A.; resources, C.A.; data curation, B.S.; writing—original draft preparation, S.B.; writing—review and editing, S.B. and K.S. and M.O. and V.S. and M.H.K and EH.C and C.A.; visualization, K.S. and M.O. and V.S. and C.A.; supervision, K.S. and M.O. and V.S. and C.A.; project administration, C.A.; funding acquisition, C.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has been financially supported by the OCP group foundation and the Mohammed VI Polytechnic University (UM6P). The lead author received financial support from the UM6P. The Horizon 2020 ACCWA project (grant agreement # 823965) in the context of Marie Skłodowska-Curie research and the innovation staff exchange (RISE) program, and the project PRIMA-S2-ALTOS-2018 "Managing water resources within Mediterranean agrosystems by Accounting for spatial structures and connectivity are also acknowledged.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The authors acknowledge all the technical support of those who helped in conducting the study. We also thank the academic editor and anonymous reviewers for accepting to review the earlier version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Abbreviations and symbols

The following abbreviations are used in this manuscript:
<table><tr><td>ELR</td><td>Environmental Lapse Rate.</td></tr><tr><td>DEM</td><td>Digital Elevation Model.</td></tr><tr><td>AWS</td><td>Automatic Weather Station.</td></tr><tr><td>MLR</td><td>Multiple Linear Regression.</td></tr><tr><td>SVR</td><td>Support Vector Regression.</td></tr><tr><td>Xgboost</td><td>Extreme Gradient Boosting.</td></tr><tr><td>Ta</td><td>Air temperature</td></tr><tr><td>Ta_5</td><td>ERA5-land's air temperature</td></tr><tr><td>Ta_st</td><td>Mesaured air temperature</td></tr><tr><td>Ta_5_ref</td><td>Reference air temperature based on ERA5-land's grid points elevation</td></tr><tr><td>Ta_5_corr</td><td>Machine learning based corrected ERA5-land's air temperature</td></tr><tr><td>ELR_cst</td><td>Constant ELR of a value of -6.5 °C/km</td></tr><tr><td>ELR_E5</td><td>Corresponding ERA5-Land ELR</td></tr><tr><td>ELR_st</td><td>Measured ELR</td></tr><tr><td>ELR_corr</td><td>Corrected ELR based on ERA5-Land corrected air temperature</td></tr><tr><td>Ta_disagg_cst</td><td>Downscaled ERA5-Land air temperaure based on constant ELR</td></tr><tr><td>Ta_disagg_MM</td><td>Downscaled ERA5-Land air temperaure based on MicorMet model</td></tr><tr><td>Ta_disagg_ML</td><td>Downscaled ERA5-Land air temperaure based on Machine learning models</td></tr><tr><td>MBE</td><td>Mean Bias Error.</td></tr><tr><td>RMSE</td><td>Root Mean Squared Error.</td></tr><tr><td>PCC</td><td>Pearson Correlation Coefficient.</td></tr><tr><td>E5</td><td>ERA5-Land grid point elevation.</td></tr></table>

554

## References

555

1. Maraun, D.; Wetterhall, F.; Ireson, A.; Chandler, R.; Kendon, E.; Widmann, M.; Brienen, S.; Rust, H.; Sauter, T.; Themeßl, M.; et al. Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user. *Reviews of geophysics* **2010**, *48*.

556

557

558

2. Maselli, F.; Pasqui, M.; Chirici, G.; Chiesi, M.; Fibbi, L.; Salvati, R.; Corona, P. Modeling primary production using a 1 km daily meteorological data set. *Climate research* **2012**, *54*, 271–285.

559

560

3. Tobin, C.; Rinaldo, A.; Schaeфli, B. Snowfall limit forecasts and hydrological modeling. *Journal of hydrometeorology* **2012**, *13*, 1507–1519.

561

562

4. Behnke, R.; Vavrus, S.; Allstadt, A.; Albright, T.; Thogmartin, W.E.; Radeloff, V.C. Evaluation of downscaled, gridded climate data for the conterminous United States. *Ecological applications* **2016**, *26*, 1338–1351.

563

564

5. Hewitt, C.D.; Stone, R.C.; Tait, A.B. Improving the use of climate information in decision-making. *Nature Climate Change* **2017**, *7*, 614–616.

565

566

6. Bjorkman, A.D.; Myers-Smith, I.H.; Elmendorf, S.C.; Normand, S.; Rüger, N.; Beck, P.S.; Blach-Overgaard, A.; Blok, D.; Cornelissen, J.H.C.; Forbes, B.C.; et al. Plant functional trait change across a warming tundra biome. *Nature* **2018**, *562*, 57–62.

567

568

7. Trisos, C.H.; Merow, C.; Pigot, A.L. The projected timing of abrupt ecological disruption from climate change. *Nature* **2020**, *580*, 496–501.

569

570

8. Fick, S.E.; Hijmans, R.J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *International journal of climatology* **2017**, *37*, 4302–4315.

571

572

9. Karger, D.N.; Conrad, O.; Böhner, J.; Kawohl, T.; Kreft, H.; Soria-Auza, R.W.; Zimmermann, N.E.; Linder, H.P.; Kessler, M. Climatologies at high resolution for the earth's land surface areas. *Scientific data* **2017**, *4*, 1–20.

573

574

10. Abatzoglou, J.T.; Dobrowski, S.Z.; Parks, S.A.; Hegewisch, K.C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Scientific data* **2018**, *5*, 1–12.

575

576

11. Navarro-Racines, C.; Tarapues, J.; Thornton, P.; Jarvis, A.; Ramirez-Villegas, J. High-resolution and bias-corrected CMIP5 projections for climate change impact assessments. *Scientific data* **2020**, *7*, 1–14.

577

578

12. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D.; et al. The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society* **2020**, *146*, 1999–2049.

579

580

13. Gelaro, R.; McCarty, W.; Suárez, M.J.; Todling, R.; Molod, A.; Takacs, L.; Randles, C.A.; Darmenov, A.; Bosilovich, M.G.; Reichle, R.; et al. The modern-era retrospective analysis for research and applications, version 2 (MERRA-2). *Journal of climate* **2017**, *30*, 5419–5454.

581

582

583

14. Saha, S.; Moorthi, S.; Wu, X.; Wang, J.; Nadiga, S.; Tripp, P.; Behringer, D.; Hou, Y.T.; Chuang, H.y.; Iredell, M.; et al. The NCEP climate forecast system version 2. *Journal of climate* **2014**, *27*, 2185–2208.

584

585

15. Kalnay, E.; Kanamitsu, M.; Kistler, R.; Collins, W.; Deaven, D.; Gandin, L.; Iredell, M.; Saha, S.; White, G.; Woollen, J.; et al. The NCEP/NCAR 40-year reanalysis project. *Bulletin of the American meteorological Society* **1996**, *77*, 437–472.

586

587

16. Kistler, R.; Kalnay, E.; Collins, W.; Saha, S.; White, G.; Woollen, J.; Chelliah, M.; Ebisuzaki, W.; Kanamitsu, M.; Kousky, V.; et al. The NCEP-NCAR 50-year reanalysis: monthly means CD-ROM and documentation. *Bulletin of the American Meteorological society* **2001**, *82*, 247–268.

588

589

590
17. Kobayashi, S.; Ota, Y.; Harada, Y.; Ebita, A.; Moriya, M.; Onoda, H.; Onogi, K.; Kamahori, H.; Kobayashi, C.; Endo, H.; et al. The JRA-55 reanalysis: general specifications and basic characteristics. *Journal of the Meteorological Society of Japan. Ser. II* **2015**, 93, 5–48.

591

592

18. Muñoz-Sabater, J.; Dutra, E.; Agustí-Panareda, A.; Albergel, C.; Arduini, G.; Balsamo, G.; Boussetta, S.; Choulga, M.; Harrigan, S.; Hersbach, H.; et al. ERA5-Land: A state-of-the-art global reanalysis dataset for land applications. *Earth System Science Data* **2021**, 13, 4349–4383.

593

594

595

19. Holden, Z.A.; Abatzoglou, J.T.; Luce, C.H.; Baggett, L.S. Empirical downscaling of daily minimum air temperature at very fine resolutions in complex terrain. *Agricultural and Forest Meteorology* **2011**, 151, 1066–1073.

596

597

20. Zhang, H.; Pu, Z.; Zhang, X. Examination of errors in near-surface temperature and wind from WRF numerical simulations in regions of complex terrain. *Weather and Forecasting* **2013**, 28, 893–914.

598

599

21. Le Roux, R.; Katurji, M.; Zawar-Reza, P.; Quénol, H.; Sturman, A. Comparison of statistical and dynamical downscaling results from the WRF model. *Environmental Modelling & Software* **2018**, 100, 67–73.

600

601

22. Alessi, M.J.; DeGaetano, A.T. A comparison of statistical and dynamical downscaling methods for short-term weather forecasts in the US N ortheast. *Meteorological Applications* **2021**, 28, e1976.

602

603

23. Zhang, G.; Zhu, S.; Zhang, N.; Zhang, G.; Xu, Y. Downscaling hourly air temperature of WRF simulations over complex topography: A case study of Chongli District in Hebei Province, China. *Journal of Geophysical Research: Atmospheres* **2022**, 127, e2021JD035542.

604

605

606

24. Vrac, M.; Drobinski, P.; Merlo, A.; Herrmann, M.; Lavaysse, C.; Li, L.; Somot, S. Dynamical and statistical downscaling of the French Mediterranean climate: uncertainty assessment. *Natural Hazards and Earth System Sciences* **2012**, 12, 2769–2784.

607

608

25. Vigaud, N.; Vrac, M.; Caballero, Y. Probabilistic downscaling of GCM scenarios over southern India. *International journal of climatology* **2013**, 33, 1248–1263.

609

610

26. Dulière, V.; Zhang, Y.; Salathé Jr, E.P. Extreme precipitation and temperature over the US Pacific Northwest: A comparison between observations, reanalysis data, and regional models. *Journal of Climate* **2011**, 24, 1950–1964.

611

612

27. Wang, J.; Fonseca, R.M.; Rutledge, K.; Martín-Torres, J.; Yu, J. A hybrid statistical-dynamical downscaling of air temperature over Scandinavia using the WRF model. *Advances in Atmospheric Sciences* **2020**, 37, 57–74.

613

614

28. Dutra, E.; Muñoz-Sabater, J.; Boussetta, S.; Komori, T.; Hirahara, S.; Balsamo, G. Environmental lapse rate for high-resolution land surface downscaling: an application to ERA5. *Earth and Space Science* **2020**, 7, e2019EA000984.

615

616

29. Ekström, M.; Grose, M.R.; Whetton, P.H. An appraisal of downscaling methods used in climate change research. *Wiley Interdisciplinary Reviews: Climate Change* **2015**, 6, 301–319.

617

618

30. Soares, P.M.; Cardoso, R.M.; Miranda, P.; de Medeiros, J.; Belo-Pereira, M.; Espirito-Santo, F. WRF high resolution dynamical downscaling of ERA-Interim for Portugal. *Climate dynamics* **2012**, 39, 2497–2522.

619

620

31. Aitken, M.L.; Kosović, B.; Mirocha, J.D.; Lundquist, J.K. Large eddy simulation of wind turbine wake dynamics in the stable boundary layer using the Weather Research and Forecasting Model. *Journal of Renewable and Sustainable Energy* **2014**, 6, 033137.

621

622

32. Laprise, R.; De Elia, R.; Caya, D.; Biner, S.; Lucas-Picher, P.; Diaconescu, E.; Leduc, M.; Alexandru, A.; Separovic, L. Challenging some tenets of regional climate modelling. *Meteorology and Atmospheric Physics* **2008**, 100, 3–22.

623

624

33. Warrach-Sagi, K.; Schwitalla, T.; Wulfmeyer, V.; Bauer, H.S. Evaluation of a climate simulation in Europe based on the WRF–NOAH model system: precipitation in Germany. *Climate Dynamics* **2013**, 41, 755–774.

625

626

34. Pan, B. Application of XGBoost algorithm in hourly PM2.5 concentration prediction. In *Proceedings of the IOP conference series: earth and environmental science*. IOP publishing, 2018, Vol. 113, p. 012127.

627

628

35. Cao, B.; Gruber, S.; Zhang, T. REDCAPP (v1.0): Parameterizing valley inversions in air temperature data downscaled from reanalyses. *Geoscientific Model Development* **2017**, 10, 2905–2923.

629

630

36. Winstral, A.; Jonas, T.; Helbig, N. Statistical downscaling of gridded wind speed data using local topography. *Journal of Hydrometeorology* **2017**, 18, 335–348.

631

632

37. Stahl, K.; Moore, R.; Floyer, J.; Asplin, M.; McKendry, I. Comparison of approaches for spatial interpolation of daily air temperature in a large region with complex topography and highly variable station density. *Agricultural and forest meteorology* **2006**, 139, 224–236.

633

634

635

38. Overland, J.E.; Wang, M. Recent extreme Arctic temperatures are due to a split polar vortex. *Journal of Climate* **2016**, 29, 5609–5616.

636

637

39. Jylhä, K.; Tuomenvirta, H.; Ruosteenoja, K. Climate change projections for Finland during the 21st century. *Boreal Environment Research* **2004**, 9, 127–152.

638

639

40. Hanssen-Bauer, I.; Achberger, C.; Benestad, R.; Chen, D.; Førland, E. Statistical downscaling of climate scenarios over Scandinavia. *Climate Research* **2005**, 29, 255–268.

640

641

41. Hofstra, N.; Haylock, M.; New, M.; Jones, P.; Frei, C. Comparison of six methods for the interpolation of daily, European climate data. *Journal of Geophysical Research: Atmospheres* **2008**, 113.

642

643

42. Tripathi, S.; Srinivas, V.; Nanjundiah, R.S. Downscaling of precipitation for climate change scenarios: a support vector machine approach. *Journal of hydrology* **2006**, 330, 621–640.

644

645

43. Pardo-Igúzquiza, E.; Chica-Olmo, M.; Atkinson, P.M. Downscaling cokriging for image sharpening. *Remote Sensing of Environment* **2006**, 102, 86–98.

646

647

44. Rodriguez-Galiano, V.; Pardo-Igúzquiza, E.; Sanchez-Castillo, M.; Chica-Olmo, M.; Chica-Rivas, M. Downscaling Landsat 7 ETM+ thermal imagery using land surface temperature and NDVI images. *International Journal of Applied Earth Observation and Geoinformation* **2012**, 18, 515–527.

648

649

649

650

651

652

653

654

655
45. Ho, H.C.; Knudby, A.; Sirovyak, P.; Xu, Y.; Hodul, M.; Henderson, S.B. Mapping maximum urban air temperature on hot summer days. *Remote Sensing of Environment* **2014**, *154*, 38–45.

650

651

46. Waichler, S.R.; Wigmosta, M.S. Development of hourly meteorological values from daily data and significance to hydrological modeling at HJ Andrews Experimental Forest. *Journal of Hydrometeorology* **2003**, *4*, 251–263.

652

653

47. Sourp, L.; Gascoin, S.; Wassim Baba, M.; Deschamps-Berger, C. Development of a snow reanalysis pipeline using downscaled ERA5 data: application to Mediterranean mountains. In Proceedings of the EGU General Assembly Conference Abstracts, 2022, pp. EGU22–5117.

654

655

656

48. Liston, G.E.; Elder, K. A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). *Journal of Hydrometeorology* **2006**, *7*, 217–234.

657

658

49. Liston, G.E.; Elder, K. A distributed snow-evolution modeling system (SnowModel). *Journal of Hydrometeorology* **2006**, *7*, 1259–1276.

659

660

50. Boudhar, A.; Duchemin, B.; Hanich, L.; Chaponnière, A.; Maisongrande, P.; Boulet, G.; Stitou, J.; Chehbouni, A. Analysis of snow cover dynamics in the Moroccan High Atlas using SPOT-VEGETATION data. *Science et changements planétaires/Sécheresse* **2007**, *18*, 278–288.

661

662

51. Chaponnière, A.; Boulet, G.; Chehbouni, A.; Aresmouk, M. Understanding hydrological processes with scarce data in a mountain environment. *Hydrological Processes: An International Journal* **2008**, *22*, 1908–1921.

663

664

52. Chehbouni, A.; Escadafal, R.; Duchemin, B.; Boulet, G.; Simonneaux, V.; Dedieu, G.; Mougenot, B.; Khabba, S.; Kharrou, H.; Maisongrande, P.; et al. An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions: The SUDMED Programme. *International Journal of Remote Sensing* **2008**, *29*, 5161–5181.

665

666

667

53. Driouech, F.; Déqué, M.; Mokssit, A. Numerical simulation of the probability distribution function of precipitation over Morocco. *Climate dynamics* **2009**, *32*, 1055–1063.

668

669

54. Bouras, E.H.; Jarlan, L.; Er-Raki, S.; Balaghi, R.; Amazirh, A.; Richard, B.; Khabba, S. Cereal yield forecasting with satellite drought-based indices, weather data and regional climate indices using machine learning in Morocco. *Remote Sensing* **2021**, *13*, 3101.

670

671

672

55. Jarlan, L.; Khabba, S.; Er-Raki, S.; Le Page, M.; Hanich, L.; Fakir, Y.; Merlin, O.; Mangiarotti, S.; Gascoin, S.; Ezzahar, J.; et al. Remote sensing of water resources in semi-arid Mediterranean areas: The joint international laboratory TREMA. *International Journal of Remote Sensing* **2015**, *36*, 4879–4917.

673

674

675

56. Dodson, R.; Marks, D. Daily air temperature interpolated at high spatial resolution over a large mountainous region. *Climate research* **1997**, *8*, 1–20.

676

677

678

57. Muñoz-Sabater, J.; Lawrence, H.; Albergel, C.; Rosnay, P.; Isaksen, L.; Mecklenburg, S.; Kerr, Y.; Drusch, M. Assimilation of SMOS brightness temperatures in the ECMWF Integrated Forecasting System. *Quarterly Journal of the Royal Meteorological Society* **2019**, *145*, 2524–2548.

679

680

58. Liu, Y.; Mu, Y.; Chen, K.; Li, Y.; Guo, J. Daily activity feature selection in smart homes based on pearson correlation coefficient. *Neural Processing Letters* **2020**, *51*, 1771–1787.

681

682

59. Sachindra, D.; Huang, F.; Barton, A.; Perera, B. Least square support vector and multi-linear regression for statistically downscaling general circulation model outputs to catchment streamflows. *International Journal of Climatology* **2013**, *33*, 1087–1106.

683

684

60. Helsel, D.R.; Hirsch, R.M. *Statistical methods in water resources*; Vol. 49, Elsevier, 1992.

685

686

61. Boser, B.E.; Guyon, I.M.; Vapnik, V.N. A training algorithm for optimal margin classifiers. In Proceedings of the Proceedings of the fifth annual workshop on Computational learning theory, 1992, pp. 144–152.

687

688

62. Cortes, C.; Vapnik, V. Support-vector networks. *Machine learning* **1995**, *20*, 273–297.

689

63. Gunn, S.R.; et al. Support vector machines for classification and regression. *ISIS technical report* **1998**, *14*, 5–16.

690

64. Vapnik, V. *The nature of statistical learning theory*; Springer science & business media, 1999.

691

65. Cristianini, N.; Shawe-Taylor, J.; et al. *An introduction to support vector machines and other kernel-based learning methods*; Cambridge university press, 2000.

692

66. Schölkopf, B.; Smola, A.J.; Bach, F.; et al. *Learning with kernels: support vector machines, regularization, optimization, and beyond*; MIT press, 2002.

693

694

67. Chen, T.; He, T.; Benesty, M.; Khotilovich, V.; Tang, Y.; Cho, H.; Chen, K.; et al. Xgboost: extreme gradient boosting. *R package version* **0.4-2** **2015**, *1*, 1–4.

695

696

68. Pesantez-Narvaez, J.; Guillen, M.; Alcañiz, M. Predicting motor insurance claims using telematics data—XGBoost versus logistic regression. *Risks* **2019**, *7*, 70.

697

698

69. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-learn: Machine learning in Python. *The Journal of machine Learning research* **2011**, *12*, 2825–2830.

699

700

70. Stone, M. Cross-validation: A review. *Statistics: A Journal of Theoretical and Applied Statistics* **1978**, *9*, 127–139.

701

71. Mohr, M.; Tveito, O. Daily temperature and precipitation maps with 1 km resolution derived from Norwegian weather observations. In Proceedings of the Proceedings 13th Conference on Mountain Meteorology/17th Conference on Applied Climatology. Citeseer, 2008, pp. 11–15.

702

703

704

72. Sluiter, R. *Interpolation methods for the climate atlas*; KNMI, 2012.

705

73. Hengl, T.; Heuvelink, G.; Perčec Tadić, M.; Pebesma, E.J. Spatio-temporal prediction of daily temperatures using time-series of MODIS LST images. *Theoretical and applied climatology* **2012**, *107*, 265–277.

706

707
74. Aalto, J.; Pirinen, P.; Heikkinen, J.; Venäläinen, A. Spatial interpolation of monthly climate data for Finland: comparing the performance of kriging and generalized additive models. *Theoretical and applied climatology* **2013**, *112*, 99–111.

75. Minder, J.R.; Mote, P.W.; Lundquist, J.D. Surface temperature lapse rates over complex terrain: Lessons from the Cascade Mountains. *Journal of Geophysical Research: Atmospheres* **2010**, *115*.

76. Shen, Y.J.; Shen, Y.; Goetz, J.; Brenning, A. Spatial-temporal variation of near-surface temperature lapse rates over the Tianshan Mountains, central Asia. *Journal of Geophysical Research: Atmospheres* **2016**, *121*, 14–006.

77. Wang, Y.; Wang, L.; Li, X.; Chen, D. Temporal and spatial changes in estimated near-surface air temperature lapse rates on Tibetan Plateau. *International Journal of Climatology* **2018**, *38*, 2907–2921.

78. Jobst, A.M.; Kingston, D.G.; Cullen, N.J.; Sirguey, P. Combining thin-plate spline interpolation with a lapse rate model to produce daily air temperature estimates in a data-sparse alpine catchment. *International Journal of Climatology* **2017**, *37*, 214–229.

79. Pepin, N.C. The possible effects of climate change on the spatial and temporal variation of the altitudinal temperature gradient and the consequences for growth potential in the uplands of northern England. PhD thesis, Durham University, 1994.

80. Shuttleworth, W.J. *Terrestrial hydrometeorology*; John Wiley & Sons, 2012.

81. Li, Y.; Zeng, Z.; Zhao, L.; Piao, S. Spatial patterns of climatological temperature lapse rate in mainland China: A multi-time scale investigation. *Journal of Geophysical Research: Atmospheres* **2015**, *120*, 2661–2675.

82. Kunkel, K.E. Simple procedures for extrapolation of humidity variables in the mountainous western United States. *Journal of Climate* **1989**, *2*, 656–669.

83. Koch, S.E.; DesJardins, M.; Kocin, P.J. An interactive Barnes objective map analysis scheme for use with satellite and conventional data. *Journal of Applied Meteorology and Climatology* **1983**, *22*, 1487–1503.

84. Kato, T. Prediction of photovoltaic power generation output and network operation. In *Integration of Distributed Energy Resources in Power Systems*; Elsevier, 2016; pp. 77–108.