Assessment of MERRA-2 and ERA5 to Model the Snow Water Equivalent in the High Atlas (1981–2019)

Abstract

Melt water runoff from seasonal snow in the High Atlas range is an essential water resource in Morocco. However, there are only few meteorological stations in the high elevation areas and therefore it is challenging to estimate the distribution of snow water equivalent (SWE) based only on in situ measurements. In this work we assessed the performance of ERA5 and MERRA-2 climate reanalysis to compute the spatial distribution of SWE in the High Atlas. We forced a distributed snowpack evolution model (SnowModel) with downscaled ERA5 and MERRA-2 data at 200 m spatial resolution. The model was run over the period 1981 to 2019 (37 water years). Model outputs were assessed using observations of river discharge, snow height and MODIS snow-covered area. The results show a good performance for both MERRA-2 and ERA5 in terms of reproducing the snowpack state for the majority of water years, with a lower bias using ERA5 forcing.

1. Introduction

The Mediterranean region is one of the areas of the world that is most impacted by water scarcity due the growing human water demand and the limited accessibility to the water resources [1,2,3]. In many countries of the Mediterranean region, meltwater runoff from the seasonal snowpack that forms in mountainous areas is of great hydrological importance [2,3,4,5]. This is particularly the case in the High Atlas Mountains in Morocco where snowmelt was estimated to contribute between 15% and 45% of annual river flow in the High Atlas catchments [6]. Snowmelt also contributes to recharge groundwater [7]. Both surface runoff and groundwater are heavily exploited to irrigate crops in the semiarid plain of Morocco [8,9].

Despite its recognized importance, key physical properties of the High Atlas snowpack are still poorly known. In particular, the main variable of interest for hydrologists and water resource planners is the snow water equivalent (SWE). Assessing the spatial distribution of SWE in any mountain region is challenging, since the snowpack is sensitive to multiple meteorological variables such as precipitation, air temperature, humidity, radiation, wind speed, etc. which are notoriously difficult to estimate in mountainous regions [10]. In addition, there is currently no spaceborne sensor which can directly measure the SWE [10]. Previous studies have been conducted in the Atlas mountains to characterize the spatial and temporal variability of the snow cover area from remote sensing [11,12,13,14]. Their contribution is important to understand the snowpack evolution, but does not allow to to have an insight into the key variables of hydrology: the distribution of SWE and snowmelt. Boudhar et al. [15], Bouamri et al. [16] computed the SWE and melt over multiple years using in situ measurements but only at the point scale.

This lack of knowledge is mainly due to the scarcity and even absence of automatic weather stations (AWS) in the mid and high elevation areas. In addition, the difficult access to some high mountain regions of Morocco makes it challenging to carry out regular measurements. To solve this issue, Baba et al. [17], Baba et al. [18], Tuel et al. [19] used meteorological reanalyses instead of weather station data to simulate the snowpack evolution over specific catchments. However, the spatial resolution of current meteorological reanalyses is too coarse to represent the influence of topography on meteorological variables, which is crucial for snow modeling in mountainous regions [20]. In order to compensate for this gap, a downscaling approach is needed to derive reanalysis datasets in a fine spatial resolution using ground observations and topographic proprieties [17,18]. Hence reanalysis data were downscaled to resolutions of 250 m [17,18] or 1 km [19]. In the High Atlas, Baba et al. [18] suggested that model resolutions coarser than 500 m may not allow a correct representation of the snowpack energy balance as the influence of the terrain slope on the radiation budget is lost below such a resolution.

ERA5 [21] and MERRA-2 [22] are the two major state-of-the-art climate reanalysis that are currently freely available at global scale. Downscaled MERRA-2 or ERA5 data were already used to feed snowpack models in various mountain regions [17,23,24,25]. Both reanalyses provide multiple meteorological variables at the hourly timestep, with a spatial resolution (MERRA-2: ∼50 km, ERA5: ∼30 km) and temporal coverage (MERRA-2: 1980 onwards, ERA5: 1950 onwards). These characteristics make them useful for long term analysis of the snow cover variability in the past decades. In addition, both reanalysis are published with a short latency of two to three weeks (MERRA-2) and 5 days (ERA5) of real time, which makes them suitable for operational planning in snow-dominated catchments as suggested by Baba et al. [26]. However, to our best knowledge, there is no study comparing the performance of both datasets over the Atlas or even over a larger region like Morocco, North Africa or the Mediterranean. This assessment is important to decide which reanalysis should be used for the implementation of a snow data assimilation scheme in Morocco [26]. It may also be useful for many other applications like the characterization of droughts [27], crop yields prediction [28] or the assessment of wind power potential [29].

In this study, we simulated the snowpack evolution in the High Atlas mountains (Tensift catchment) from 1981 to 2019 using a distributed snowpack evolution model [30] [SnowModel] forced by ERA5 and MERRA-2 meteorological data. The results were evaluated using available in situ meteorological data, MODIS-derived snow-covered fraction (SCF) and observed discharge. We also evaluate SWE directly from the recently released ERA5-Land. ERA5-Land is a new land reanalysis dataset providing an estimation of several land variables at 0.1° × 0.1° of spatial resolution, and hourly temporal resolution [31].

2. Study Area and Data Source

2.1. Study Area

The study area comprises the upper catchments of the Tensift river basin on the northern slope of the High Atlas, including the Toubkal peak (4167 m.a.s.l), the highest summit in North Africa (Figure 1). The seasonal snow covers the study area in the middle and high altitudes above 1500 m.a.s.l in winter [6,13,18]. The studied sub-catchments of the Tensift basin provide a major source of freshwater for the Haouz plain, where the average annual precipitation varies between 150 mm and 200 mm, while it reaches 800 mm in high altitudes regions [6,32]. While the simulation domain covers all these sub-catchments in Figure 1, the evaluation was limited to some sub-catchments depending on the availability of observation data.

2.2. Merra-2 Reanalysis

The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), is provided by the National Aeronautics and Space Administration (NASA) Global Modeling and Assimilation Office (GMAO) from 1980 onward. MERRA-2 was the first long-term global reanalysis to assimilate space-based observations of aerosols and represent their interactions with other physical processes in the climate system [33]. All of MERRA-2 atmospheric variables are provided at 0.5° × 0.625° spatial resolution with hourly time-step.

2.3. Era5 Reanalysis

ERA5 is the fifth generation of atmospheric reanalyses of the global climate produced by the European Center for Medium-Range Weather Forecasts (ECMWF) [21]. It is a key element of the EU-funded Copernicus Climate Change Service (C3S). ERA5 provides a large number of atmospheric, land and oceanic climate variables at hourly step. ERA-5 have a high spatial resolution as well, with all of its atmospheric variables are interpolated at 0.25° × 0.25° spatial resolution.

2.4. Dem

A DEM of 200 m of spatial resolution was used for this study specifically. It was resampled by using bicubic interpolation from the Shuttle Radar Topography Mission Void Filled (SRTM-Void Filled Version 4.1) [34]. The choice of this resolution was based on a recent research in the High Atlas that suggest that 200–250 m resolutions are the best trade-off between the computational time and the reliability of the snowpack evolution simulation in the High Atlas [18].

2.5. Land Cover

The land cover map was derived from the European Space Agency Climate Change Initiative (ESA-CCI) classification map [35]. The bare soil is the predominant class over the study area.

2.6. Modis Snow Cover Area

We used an 18-year time series of daily gap-filled snow cover fraction (SCF) maps that were derived from MODIS snow products (NASA MOD10A1 collection 6). Each map provides the snow cover fraction per pixel with a spatial resolution of 500 m [36]. The method to retrieve SCF maps from MODIS and to interpolate cloud-masked pixels was described in [13,17]. We have computed the snow-covered area for each sub-catchment by using SCF products.

2.7. Hydro-Climatic Observations

We used annual river discharge data at the outlets of five rivers (

Table 1. The outlet of each catchment.

Catchment	Rdat	Zat	Ourika	Rheraya	Nfis
Outlet	Sidi Rahal	Taferiat	Aghbalou	Tahanaout	Iguir Nkouris
Area (Km${}^2)$	552	525.9	502.6	225.3	823.8
Max elevation (m)	3476	3847	3996	4098	4088
Mean elevation (m)	1718.7	1831.7	2444.5	2165.8	2202.1

). The discharge observations are available from 2000 to 2015.

2.8. Era5-Land Swe

We used daily SWE reanalysis datasets generated by the European Centre for Medium-Range Weather Forecasts (ECMWF) [31]. The SWE reanalysis are available from 1981 to present at 0.1° × 0.1° spatial resolution.

3. Methodology

3.1. Snowpack Evolution Modeling

We used the spatially distributed snow evolution modeling system (SnowModel) [30]. SnowModel is composed of 4 sub-models:

• Micromet computes the spatial distribution of meteorological data (precipitation, air temperature, relative humidity, wind speed, incoming solar radiation and incoming longwave radiation) over the model grid. This distribution can be done at daily, 3 hourly or daily scale depending on the temporal resolution of the input [37].
• EnBal computes the energy fluxes between the snowpack and the atmosphere. It simulates the surface temperature, internal energy, net radiation, sensible heat, latent heat, etc., using the MicroMet outputs. In addition, it simulates the latent and sensible heat fluxes and determines the potential amount of melt [30].
• Snowpack computes the evolution of the height of the snowpack (HS) based on the MicroMet and EnBal outputs.
• SnowTran-3D, this sub-model computes the redistribution of the snow due to wind transport and the blowing snow sublimation [38].

SnowModel has previously been tested with success in a wide range of climates [39], including semi-arid and Mediterranean mountains [17,40,41]. In this study we used the version: SnowModel_2020_05_14.

To spatialize the meteorological data from ERA5 and MERRA-2 centroids, a horizontal interpolation is performed using the Barnes objective analysis scheme, which applies a Gaussian distance-dependent weighting function to interpolate the station data to the model grid resolution. Thereafter, a horizontal distribution is applied to take into account the effect of elevation, slope and aspect on the gridded meteorological forcing. In the case of air temperature, a linear lapse rate correction is applied, while for the precipitation a non-linear lapse rate is applied for the correction. For further description of this distribution the reader can be referred to [17,30].

3.2. Evaluation Methodology

The first evaluation is based on a comparison of snow cover simulated by SnowModel and the MODIS snow cover data [42,43]. The simulated SWE at every day (MODIS revisit time) was converted to binary snow cover area (SCA) using a SWE${}_0$ threshold. This means that if the SWE value of a cell is greater than SWE${}_0$, then the SCA is equal to 1, otherwise the SCA is equal to 0. SWE${}_0$ is optimized by Baba et al. [17] and set to 20 mm. The simulated SCA maps were converted to SCF by resampling to MODIS SCF resolution (Figure 2). The resampling method is the mean resampling method. For MODIS it provides snow cover fraction (SCF) for each pixel.

Then, we computed the daily mean SCF over the entire study area and compared it to the MODIS daily main SCF. This examination was conducted at the scale of each sub-catchment from 2000 to 2018. In addition, Snow cover duration which represents the number of days with snow cover is computed from both MODIS SCF and modeled SCF.

We also evaluated the simulations using the daily HS records observed at Oukaimeden station (3200 m.a.s.l) during 7 snow seasons. The evaluation was carried out by comparing the variations of each HS and by computing the root mean square error (RMSE) between modeled and observed HS.

SnowModel outputs were further assessed by comparing the total surface water input (defined as the liquid water leaving the snowpack, equal to rainfall rate in snow-free areas) over each sub-catchment to the observed discharge. We performed this comparison at a yearly time step and qualitatively only since several hydrological processes which are involved to convert the surface water input to discharge are not taken into account in this study. The comparison consists on comparing the variation of the observed discharge with 30–50% of the modeled surface water input (the sum of the precipitation and snowmelt that could contributes to the discharge).

Moreover, we evaluated the modeled precipitation at Oukaimeden based on precipitation observation, and the air temperature at four different stations (CAF, Neltner, Oukaimeden and Tahanaout). They are all located at Rheraya catchment. Ground observations are from a research weather station that is not part of an operational network hence are not assimilated in ERA5 or MERRA-2 reanalyses.

Finally, we have compared the daily SWE derived from ERA5-Land to the modeled SWE. The comparison was performed by comparing the evolution of each SWE by computing the root mean square error (RMSE) between modeled and observed total SWE (mm) in the study area. The main purpose of this comparison is to assess the performance of ERA5-Land in terms of SWE.

These methods of evaluation are summarized in Figure 2.

4. Results

4.1. Hs Evaluation

The evolution of modeled and observed daily snow height at Oukaimeden is shown in Figure 3. The evaluation of this evolution generally shows that SnowModel forced by both ERA5 and MERRA-2 performs well in terms of event detection (snowfall and melt). To complete the comparison we present the differences between modeled and observed HS in

Table 2. Comparisons between modeled and observed HS (cm) (RMSE and mean error).

RMSE	2009	2010	2013	2014	2015	2016	2017
ERA5	24.46	43.71	32.01	21.83	15.10	28.94	26.70
MERRA-2	25.53	26.29	40.46	28.34	24.79	26.99	34.93
Mean Error	2009	2010	2013	2014	2015	2016	2017
ERA5	−11.35	38.58	14.09	4.62	1.25	2.57	5.45
MERRA-2	18.64	4.46	−28.74	6.52	12.65	8.91	18.69

. The mean RMSE compared to observed HS is equal to 27.76 cm for ERA5 simulations and 29.61 cm for MERRA-2 simulations. In some seasons (e.g., 2010 and 2017) the model did not perform well. For these years, the RMSE is respectively.equal to 43.71 cm, 26.70 cm (ERA5), 26.29 cm, 34.93 cm (MERRA-2). Therefore in most cases ERA5 performs better than MERRA-2 in modeling the snow height evolution at Oukaimeden.

4.2. Scf Evaluation

Considering the SCF evolution with respect to MODIS data (Figure 4), we find that SnowModel driven by ERA5 outperforms SnowModel forced by MERRA-2. This is confirmed by statistical metrics, such as the the adjusted R-squared (R${}^2$) of the mean of daily SCF R${}^2$ (MODIS, ERA5) and R${}^2$ (MODIS, MERRA-2) are respectively equal to 0.82 and 0.66. Moreover, the same indices calculated at catchment scale show a similarly good performance of both simulation with a slight superiority of ERA5 simulation (

Table 3. Mean R${}^2$ scores from 2000 to 2018 between modeled and observed SCF for each sub-catchment area of Tensift.

Catchment	Ghdat	Nfis	Ourika	Rheraya	Zat	Seksawa	Mhand	Assif
Mean R${}^2$ ERA5	0.83	0.77	0.82	0.81	0.74	0.74	0.79	0.89
Mean R${}^2$ MERRA-2	0.70	0.67	0.85	0.79	0.65	0.67	0.75	0.76

).

The comparison of the snow cover duration (SCD) from the model output and the MODIS data shows that both model configurations reproduce well the duration of snow cover with a mean error of 15% (ERA5) and 10% (MERRA-2). It is in line with the mean error between modeled HS and observed HS at Oukaimeden (Table 2), the bias is generally positive.

4.3. Discharge Evaluation

The comparison of the observed annual mean instantaneous discharge (m${}^3$/s) with 30–50% of modeled surface water input shows that the models reproduce well the inter-annual variability of the annual surface water input with a correlation varying from 0.73 to 0.87 for ERA5 simulations and 0.65 to 0.82 for MERRA-2 (Figure 5 and Table 3) Some water years are poorly simulated (e.g., WY 2009 for NFIS, WY 2012 for Rheraya (Figure 5)). We also observe that the magnitude of the inter-annual variability of the annual surface water input is quite biased.

4.4. In Situ Precipitation

The comparison of the evolution of cumulative precipitation at Oukaimeden shows that the RMSE between both ERA5, MERRA-2 and the observed cumulative precipitation is suitable (90 mm for ERA5 and 106 mm for MERRA-2). However, in some years such as WY 2012 a few anomalies of precipitation are and many events are not detected as shown in (Figure 6). More, we computed the Mean R${}^2$ scores between the cumulative precipitations and it is equal to 0.96 and 0.92 for ERA5 and MERRA-2 respectively.

4.5. Air Temperature

The air-temperature (T) comparison was established in four regions at the Rheraya catchment. The results confirm that the ERA5 daily average temperature (from 2002 to 2016) is closer to the observations compared to MERRA-2 (

Table 4. Indicators to evaluate MERRA-2 and ERA5 reanalysis with observed temperature (RMSE and R${}^2$).

	Oukaimeden	Tachedirt	Neltner	CAF
ERA5 RMSE (C)	1.84	2.45	2.37	2.02
MERRA-2 RMSE (C)	3.73	2.77	3.73	3.92
ERA-2 R${}^2$	0.90	0.95	0.93	0.94
MERRA-2	0.78	0.82	0.77	0.79

). These results could be confirmed visually in Figure 7. The distribution of daily observed temperature versus modeled ones for four different stations shows that the standard deviation for between observed temperature and ERA5 temperature is much lower than the standard deviation between observed temperature and MERRA-2 temperature.

4.6. Era5-Land Swe

The comparison of the evolution of total SWE at the study area shows that the mean R${}^2$ between modeled SWE and ERA5-Land SWE is equal to 0.56 for ERA5 and 0.46 for MERRA-2. Figure 8 shows the interannual variation of the annual mean SWE averaged for the High Atlas mountains from the different outputs. The ERA5-Land SWE is considerably overestimating the SWE compared with the simulations. The mean of SWE averaged for the different products is equal to 3.24 mm, 2.47 mm and 10.40 mm for modeled SWE with ERA5, modeled SWE with MERRA-2 and ERA-Land SWE respectively.

5. Discussion

This study provides an evaluation of ERA5 and MERRA-2 to force a model of the snowpack evolution in a the High Atlas of Morocco. In this study, the snowpack evolution model SnowModel was forced by a the reanalysis for the 1981–2018 period. Given the initial coarse resolution of the meteorological reanalysis, the evaluation of the different outputs shows that SnowModel performs well with ERA5 and MERRA-2. Indeed, we assume that ERA5 has a better estimation of the key meteorological data (precipitation and temperature) when using MicroMet as downscaling scheme. In this case, it subsequently impacts the estimation of the key variables of hydrology (snow water equivalent, snow depth and surface water input). Thus, even the evaluation of the snow depth and the surface water input shows a higher quality for ERA5 simulations. These results were expected since the spatial resolution of ERA5 is lesser than MERRA-2 and therefore should give a better representation and distribution of the meteorological data. It is highly encouraged to assess this hypothesis with other downscaling techniques.

Despite the good quality of the simulations with ERA5, some discrepancies were found for some years. As an example, comparing the modeled and the observed snow height at Oukaimeden in January 2009 and February 2013 (Figure 3) show an important bias. This is mainly due to the absence of certain of snowfall events. These failures of the model are normal and also expected since it is difficult to catch local events. To this date, the estimation of precipitation and snowfall by numerical models is very challenging in mountainous regions [25,44,45]. We note also that the comparison of precipitation has been carried out over very complex terrain (Oukaimeden) where the altitude varies from 2400 m to 3200 m in less than 3 km. Thus, we believe that the focus of future work should be on the correction of precipitation using data assimilation techniques [26,46].

This study shows also that the modeled snow-covered area reproduces that observed by MODIS relatively well, but with a few differences. These differences enable us to better understand some possible sources of model and forcing errors. In the January 2009 we notice that the modeled SCF is approximately equal to the half of the observed one Figure 4. Similar underestimation is shown in the HS evolution, the observed HS is approximately equal to the double simulated Figure 3 for both ERA5 and MERRA-2. Similar underestimation could be observed also in February 2013. We also note that for many years where the SCF is adequately reproduced, the snow height is also faithfully simulated at Oukaimeden. In the High Atlas there is not much vegetation and forests especially in the middle and high altitudes where the snowpack is present for several months. Therefore, we think that MODIS limitations because of the vegetation will not affect the quality of the derived SCF [47,48] . Thus we think that the model reproduce well the SCF since the validation data are reliable.

Both modeled and observed outputs highlight the large interannual variability of the snowpack evolution, in term of amount of snow and of duration of the snowpack. There are several water years that the SCF exceeds 25% km${}^2$ for more than 30 days (e.g., 2005, 2011, 2014, 2016), and some water years where the SCF barely exceeds 25% (e.g., 2009, 2012, 2018) (Figure 4. HS evolution show that in the high elevated areas (Oukaimeden 3200 m.a.s.l) the snowpack generally is evaluating between December and May of each water year, with a height varying from 50 cm to 150 cm (Figure 3). With an exception for 2015 W.Y where the snow season started in March and lasted only two month (Figure 3). This anomaly of the snow season is reproduced by both simulations of HS and SCA, and by the observations of HS and SCA. This is confirmation that the reanalysis can represent extreme events well. Therefore, it could be a good resource to better understand the snowpack evolution given the lack of in situ measurements.

We conducted an intercomparison of three SWE products: ERA5-Land SWE and modeled SWE using SnowModel with ERA5 and MERRA-2. We have noticed that the modeled SWE are relatively close in terms of the interanual variation of the mean SWE averaged for the zone of interest. While, ERA5-Land SWE is overestimating the SWE compared to the modeling SWE. Given the previously proven quality of the simulations and the high value of the annual mean SWE derived from ERA5-Land (10.40 mm) we assume that ERA-Land is overestimating the SWE in the High Atlas mountains.

6. Conclusions

We presented a methodology to assess the ability of MERRA-2 and ERA5 to be downscaled and used as forcing in an independent (offline) snowpack model (SnowModel) in the High Atlas from 1981 to 2019. The results showed both reanalysis using SnowModel allow the reproduction of the snowpack state with the ERA5-driven simulations being superior. This might be related to the fact that ERA5 data are approximately four times denser than MERRA-2 data. These results suggest that ERA5 reanalysis could be used to simulate SWE distribution in the High Atlas. This leads us to think that a good estimation of SCF has a positive impact on the estimation of HS, SWE and the other variables.

However, there remain biases in the ERA5 meteorological reanalyses, which could be significantly reduced by assimilating remote sensing snow cover area products (e.g., Sentinel-2 or MODIS) [25,49]. Indeed, assimilation snow cover area derived from remote sensing products in hydrological and snowpack evolution models has shown a great added value in the literature. Andreadis and Lettenmaier [50] updated the SWE state based on MODIS observation and they noticed that their assimilation scheme improvement was evident in lower and middle elevations, especially during melt period. Thus we believe that the assimilation will lead to a better estimate of the state of the snowpack in the High Atlas mountains. Remote sensing snow cover area products could also be used to reduce uncertainty in hydrological modelling [43,51] and on correct the estimation of snowfall spatial distribution in mountainous area [52].

At the same time, the results should be complemented with more in-situ observations. Indeed, planning two or three drone acquisition in pilot catchment (e.g Rheraya or Ourika) could provide an accurate estimation of the snow height with a distributed manner (<1 m of spatial resolution). The obtained HS maps would provide a stronger evaluation criteria of the models performance.

Finally, we compared the simulated SWE with reanalyses against ERA5-Land SWE. The results obtained show important discrepancies between the simulated outputs and ERA-Land SWE dataset. This leads us to suppose that ERA5-Land SWE is not suitable for the High Atlas mountains, and a move to more detailed modeling is recommended to better reproduce the snowpack state.