2.1 Initialization

ERA5-Land is not produced as a single continuous simulation for the entire period. The production is conducted in three independent streams, as shown in Fig. 1. To avoid or minimize discontinuities between streams, a careful initialization procedure is needed for each of them. Particular attention must be given to variables carrying long memory. As an example, Fig. 3 shows time series of deep soil moisture in a band of latitude between 60 and 20^∘ S, where the averaged annual soil moisture variability is low. This period includes several production streams of ERA5. ERA5 initialized each stream with ERA-Interim soil moisture initial conditions, which has a different climatology than ERA5. In ERA5, a 1-year spin-up was used for each production stream and it is normally long enough for atmospheric variables. However, it is not sufficient for deep soil moisture to reach equilibrium, which leads to discontinuities between two production streams, as shown in Fig. 3.

Figure 3(a) Mean differences between ERA5-Land and ERA5 soil moisture (sm) time series for the fourth soil layer of the Carbon Hydrology-Tiled ECMWF Scheme for Surface Exchanges over Land (CHTESSEL) land surface model (100–289 cm), averaged for the latitudes 60–20^∘ S. The temporal resolution is 12 h. Panel (b) is the same as (a) but showing the raw time series of ERA5 (dashed light grey curve) and ERA5-Land (dark curve).

Download

The strategy followed to initialize the ERA5-Land production stream starting in 2001 (stream-1) was to use the latest year of a long, prior ERA5 stream, and letting 3 further spin-up years to allow a long spin-up period (see Fig. 1). While this strategy provides satisfactory results for most continental masses around the world, discontinuities are still possible at areas with very low variability of soil moisture (deserts and polar regions). Particular attention was given to the treatment of permanent snow-covered regions. The current model formulation (as in ERA5) does not have an independent treatment of glaciers. Grid points with glaciers are assigned with a constant snow mass of 10 m. ERA5-Land streams are initialized on the 1 January, and a glacier mask is applied to snow mass to guarantee the correct spatial representation of glaciers. A threshold of 50 % of a grid box covered by ice is used, below which the snow depth keeps the value computed by the snow scheme of the land model. Values above the threshold assign a snow water equivalent value of 10 m. This condition is used to avoid grid points near glaciers with large unrealistic snow depth that result from the interpolation from ERA5 fields to ERA5-Land. For the stream starting in 1981 (stream-2), a long, prior ERA5 stream was not available. The strategy in this case was to initialize using a ERA5-Land climatology of 1 January for the period 2001–2018 and then allow the system to spin up for 3 years. A similar strategy was used to initialize the stream starting in 1950 (stream-3), but in this case a 1981–2010 climatology was used and only 1 spin-up year was feasible. The latter was limited by the availability of forcing data.

2.2 Static and climatological fields

As in ERA5, the land characteristics are described using several time-invariant fields. These consist of the land–sea mask, the lake cover and depth, the soil and vegetation type, and the vegetation cover. In addition, surface albedo and leaf area index are prescribed as monthly climatologies. The complete list of time invariant fields with their information source is provided in Table A1 of Appendix A.

2.3 Atmospheric forcing

ERA5-Land is driven by atmospheric forcing derived from ERA5 near-surface meteorology state and flux fields. The meteorological state fields are obtained from the lowest ERA5 model level (level 137), which is 10 m above the surface, and include air temperature, specific humidity, wind speed, and surface pressure. The surface fluxes include downward shortwave and longwave radiation and liquid and solid total precipitation. These fields are interpolated from the ERA5 resolution of about 31 km to ERA5-Land resolution of about 9 km via a linear interpolation method based on a triangular mesh. The atmospheric forcing built for ERA5-Land is hourly and consistent over the entire production period, and it is the result of the assimilation of a large number of conventional meteorological and satellite observations through a four-dimensional variational assimilation system (4D-Var) and simplified extended Kalman filter (SEKF) systems as described in Hersbach et al. (2020). Previous land reanalyses have included corrections to the precipitation forcing to address limitations of the precipitation fields of the atmospheric reanalysis. This is not the case in ERA5-Land mainly due to the (1) enhanced quality of ERA5 precipitation when compared with previous atmospheric reanalyses (e.g. Beck et al., 2019; Tarek et al., 2020; Nogueira, 2020) and (2) reduced dependencies on external data that would limit the near-real-time data availability. However, air temperature, humidity, and pressure are corrected for the altitude differences between ERA5 and ERA5-Land grids. This correction involves four steps: (i) relative humidity is computed from interpolated, uncorrected fields; (ii) air temperature is adjusted for the altitude differences using a daily environmental lapse rate (ELR) field derived from ERA5 lower troposphere temperature vertical profiles (Dutra et al., 2020); (iii) surface pressure is corrected for the altitude differences and correction of temperature; and (iv) specific humidity is computed using the corrected temperature and pressure assuming that there is no change in relative humidity. Dutra et al. (2020) present a detailed evaluation of this methodology comparing the use of a constant (time and space) ELR with daily ELR fields derived from ERA5. This methodology was shown to reduce the mean absolute error (MAE) of daily maximum temperature by 10 % and by 4 % for daily minimum temperature with respect to ERA5 when compared with 2941 stations over the western US. The importance of the ELR correction is shown in the orography map of the Alpine region around Switzerland in Fig. 4. ERA5 misses many of the highest Alpine peaks due to the coarser resolution. Taking into account these orographic differences is important for other land variables such as surface temperature. For instance, the middle and bottom rows show a more realistic spatial pattern of surface temperature in ERA5-Land, with a clear cold signal over the higher peaks. It is also remarkable to see how ERA5-Land is able to resolve lakes such as Geneva (Léman), Neuchâtel, and Constance. This is especially visible at 06:00 UTC, when the lake surface temperature is still significantly warmer than the land (Fig. 4f).

Figure 4ERA5 (a, c, e) and ERA5-Land (b, d, f) orography (a, b) of the Alpine region around Switzerland. Panels (c, d) and (e, f) show the ERA5 and ERA5-Land land surface temperature (LST) estimates on 15 July 2018 at 06:00 UTC (c, d) and 15:00 UTC (e, f), respectively. The locations of lakes Geneva (Léman), Neuchâtel, and Constance are indicated in panel (b).

2.4 Land surface model

The core of ERA5-Land is the ECMWF land surface model: the Carbon Hydrology-Tiled ECMWF Scheme for Surface Exchanges over Land (CHTESSEL). The main updates with respect to the land component of ERA-Interim are (a) a revised soil hydrology, introducing an improved formulation of the soil hydrologic conductivity and diffusivity (that now is variable as a function of soil texture) and surface runoff based on variable infiltration capacity (Balsamo et al., 2010); (b) a fully revised parametrization of the snow scheme, changing the hydrological and radiative properties of the snowpack (Dutra et al., 2010); (c) the introduction of a climatological seasonality of vegetation, in contrast to the fixed vegetation in ERA-Interim (Boussetta et al., 2013b); (d) a new scheme for bare soil evaporation, allowing soil moisture to reach values below the wilting point (Albergel et al., 2012); (e) introduction of a lake model to represent the thermodynamics of inland water bodies (Balsamo et al., 2012); and (f) a parametrization that allows the estimation of land carbon fluxes – net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_eco) – in a modular way with the Jarvis approach that computes the stomatal conductance without affecting the transpiration components (Boussetta et al., 2013a).

The land surface model version used in ERA5-Land was operational at ECMWF in 2018 with the model cycle Cy45r1. A detailed description of the model can be found in chapter IV of the Integrated Forecasting System (IFS) documentation (https://www.ecmwf.int/node/18714, last access: December 2020). Compared with the model version of ERA5, the differences are mostly technical, with the exception of (i) an updated parametrization of the soil thermal conductivity following Peters-Lidard et al. (1998) that takes into account the ice component in the case of frozen soil, (ii) a fix to improve conservation for the soil water balance, and (iii) rain over snow is accounted for and it does not accumulate in the snowpack. Potential evapotranspiration (PET) flux in ERA5 suffers from a bug present in IFS cycle 41r2 that affects PET computation over forests and deserts. This problem has been corrected in ERA5-Land, and unlike in ERA5, ERA5-Land includes PET in the portfolio of products. Given the importance of this variable for some applications, it is worth clarifying that PET is computed by making a second call to the surface energy balance assuming a vegetation type of crops and no soil moisture stress. In other words, evaporation is computed for agricultural land as if it is well watered and assuming that the atmosphere is not affected by this artificial surface condition. The latter may not always be realistic. Therefore, despite the fact that PET is meant to provide an estimate of irrigation requirements, one has to be cautious especially for arid conditions, since the method can give unrealistic results due to too-strong evaporation forced by dry air.

3.1 ERA-Interim

ERA-Interim (Dee et al., 2011) is the former ECMWF reanalysis providing estimates for the atmosphere, ocean waves, and land surface. It had supported scientific progress during the previous decade and is still widely used today by the scientific community. It includes information on multiple land and atmospheric variables, and it is available from January 1979 to August 2019. It used the IFS version Cy31r2 (more detailed information available at https://www.ecmwf.int/en/publications/ifs-documentation, last access: 27 August 2021), corresponding to the IFS 2006 release, with a spatial horizontal resolution of about 80 km and 60 levels in the vertical from the surface to 0.1 hPa. The system includes a 4D-Var system (Rabier et al., 2000) providing analyses fields at a temporal resolution of 6 and 3 h for short-forecast fields (like for precipitation and fluxes). The main difference between the land component of ERA-Interim and those of both ERA5 and ERA5-Land is that the former is based on TESSEL (van den Hurk et al., 2000), which is considered the precursor of the current CHTESSEL scheme. In ERA-Interim, the soil moisture and soil temperature analyses are based on a local optimal interpolation scheme (Mahfouf, 1991; Douville et al., 2000) that assimilates surface synoptic (SYNOP) observations of temperature and relative humidity at screen level (2 m). The snow depth analysis is independent of the soil wetness and is based on a Cressman analysis that assimilates SYNOP snow reports and snow-free satellite observations (Drusch et al., 2004).

3.2 ERA5

ERA5 is the latest comprehensive ECMWF reanalysis and has replaced ERA-Interim. It is based on a version of the ECMWF IFS (Cy41r2) that was operational in 2016. ERA5 provides hourly estimates of the global atmosphere, land surface, and ocean waves from 1950 and is updated daily with a latency of 5 d. Its state estimates are based on a high-resolution (HRES) component at a horizontal resolution of 31 km and with 137 levels in the vertical spanning from the surface up to 0.01 hPa. Information on uncertainties in these are provided by a 10-member ensemble of data assimilations (EDA) at half the horizontal resolution. Both the HRES and EDA ERA5 data assimilation use background-error estimates that utilize the output from the EDA. The land component of ERA5 is, like ERA5-Land, based on the CHTESSEL model, though at a resolution of 31 km rather than 9 km. ERA5 uses new analyses of sea-surface temperature and sea-ice concentration, variations in radiative forcing derived from CMIP5 specifications, and various new and reprocessed observational data records.

The data assimilation system consists of an incremental 4D-Var component (Courtier et al., 1994) for upper-air and near-surface components, an ocean-wave optimal interpolation scheme, and a dedicated LDAS. The LDAS comprises a two-dimensional optimal interpolation scheme for the analysis of screen-level 2 m temperature and relative humidity, and for snow (depth and density), a point-wise simplified extended Kalman filter (de Rosnay et al., 2013) for three soil moisture layers in the top 1 m of soil, and a one-dimensional optimal interpolation for soil, ice, and snow temperature. Details of the ERA5 configuration are given in Hersbach et al. (2020), which also contains a basic evaluation of characteristics and performance for the segment from 1979 onward. The performance of the component from 1950 to 1978, the back extension which was made available later, is described in Bell et al. (2021), and a detailed analysis for surface temperature and humidity is provided by Simmons et al. (2021). Technical details are also provided in the online documentation (https://confluence.ecmwf.int/display/CKB/ERA5:+data+documentation, last access: 23 August 2021).

Table 1Overview of the main characteristics of ERA-Interim, ERA5, and ERA5-Land.

^* Availability at the time of submitting this paper. n/a – not applicable.

Download Print Version | Download XLSX

Both ERA5 and ERA5-Land are produced as part of the Copernicus Climate Change Service that ECMWF operates on behalf of the European Commission and are available from the C3S Climate Data Store (CDS). ERA5 and ERA5-Land have a large and diverse user base (more than 40 000 users at the end of 2020). Table 1 summarizes the main characteristics of ERA-Interim, ERA5, and ERA5-Land.

3.3 Soil moisture

To evaluate the quality of the soil moisture estimates from the various reanalyses, a large number (>800) of in situ sensors in the period 2010 to 2018, many providing hourly measurements, were used. These sensors belong to the networks that are listed in Table 2. The networks are located in North America, Europe, Africa, and Australia, and all the observations were retrieved from the International Soil Moisture Network (Dorigo et al., 2011, 2021). Three reanalysis soil layer estimates were compared to measurements by sensors at three different depths. ERA5 and ERA5-Land top layer soil moisture estimates (0–7 cm) were compared to in situ sensors at 5 cm depth in North America, Africa, Europe, and Australia. A more in-depth study was performed in North America, where most of the sensors are located. For this region, surface soil moisture from ERA-Interim was also evaluated. In addition, ERA5, ERA5-Land, and ERA-Interim soil moisture for the second (7–28 cm) and third layers (28–100 cm) were evaluated against in situ measurements at 20 and 50 cm depth, respectively. In situ measurements were compared to the closest grid point of the ERA5, ERA5-Land, and ERA-Interim grids if the closest grid point was not farther than the respective model resolution. In situ time samples were selected in a ±1 h window with respect to the reanalysis timestamp. If several observations were retrieved in a single window, the average was computed. In order to consider an observed time series suited for the comparison, a minimum of 150 samples was required for the study period. To remove the seasonal cycle, anomaly time series were also computed. The soil moisture anomaly values at time t (SM_AN(t)) were computed from the original time series (SM(t)) computing the mean (SM‾) and the standard deviation (σ_SM) of soil moisture a ±17 d window as follows:

The following metrics were computed between the reanalyses estimates and the in situ time series: the standard deviation of the difference (STDD, equivalent to the unbiased RMSD), the bias, and the Pearson correlation coefficient (R). The latter was also computed for the anomaly time series (R_AN). The results were grouped per continent and their distribution presented in box plots. Values are considered outliers if they are greater than q75+1.5×(q75-q25) or less than q25-1.5×(q75-q25), with q₂₅ and q₇₅ the 25th and 75th percentiles, respectively.

Table 2In situ measurement networks used to evaluate soil moisture. The different columns contain the region and the name of the network, the depths of the probes used, and the bibliographic reference to the network.

Download Print Version | Download XLSX

3.4 Snow

Snow depth estimates from reanalyses were compared to two sets of observational data. The first one comprises 10 sites distributed among North America, Europe, and Japan. They were selected as reference sites to evaluate cold processes by models participating in the Earth System Model – Snow Model Intercomparison Project (ESM-SnowMIP) (Krinner et al., 2018; Ménard et al., 2019). These sites provide benchmarking data for cold processes in maritime, Alpine, and taiga types of snow cover and on different types of climates. Table 3 presents these stations.

Table 3List of ESM-SnowMIP sites used for the evaluation of the snow parameters, adapted from Krinner et al. (2018).

Download Print Version | Download XLSX

The second dataset is retrieved from the Global Historical Climatology Network-Daily (GHCN-daily) (Menne et al., 2012b) from 1 July 2010 to 30 June 2018. The version used is v3.24 (Menne et al., 2012a). This network integrates thousands of land surface stations across the globe. The daily observed snow depth product was compared to the daily averaged snow depth estimates from ERA-Interim, ERA5, and ERA5-Land at 00:00 and 12:00 UTC. To compare the reanalysis data with the snow depth observations, the prognostic snow water equivalent (SWE, in units of m water equivalent) and snow density (ρ_snow, kg m⁻³) from the reanalysis were combined to compute the actual snow depth (SD, m):

where ρ_water=1000 kg m⁻³ is the reference density of the water. The GHCN observations were only considered if snow depth values were positive and missing values were lower than 50 % of the recorded time series. Also stations where the snow depth reported is lower than 1 cm in more than 5 % of the total number of days were removed. Finally, stations located in coastal areas with more than 50 % water in the pixel and on permanent snow area (glaciers) were removed. More than 6000 stations passed the quality filters, and their locations can be found in Fig. 9. The quality of the snow depth from the reanalyses was evaluated at the hemispheric scale by computing the mean bias (here defined as reanalysis estimate minus in situ observation) and root mean square error (RMSE) for the months between December and June of the 2010–2018 period.

3.5 Lakes

In 2015, the CHTESSEL land surface scheme of the operational IFS introduced a lake tile, which represents lakes, reservoirs, rivers, and coastal (subgrid) waters, and is based on the FLake (Fresh-water Lake) model of Mironov et al. (2010a). FLake is a one-dimensional model, which uses an assumed shape for the lake temperature profile including the mixed layer (uniform distribution of temperature) and the thermocline (its upper boundary located at the mixed-layer bottom, and the lower boundary at the lake bottom). To run FLake, the lake location (or fractional cover), lake depth (most important parameter, preferably bathymetry), and lake initial conditions are required. For the best performance, lake depth should be updated with the latest available information to ensure that depths are close to observed values, as overestimated depths can be blamed for cold biases in summer temperatures or lack of ice. The state of lakes in FLake is described by seven prognostic variables: mixed-layer temperature, mixed-layer depth, bottom temperature, mean temperature of the water column, shape factor (with respect to the temperature profile in the thermocline), temperature at the ice upper surface, and ice thickness. In this paper, the lake surface water temperature (LSWT) estimates from the FLake model embedded in ERA5 and ERA5-Land were compared to in situ observations from three different sources of data, during ice-free periods:

• The Alqueva reservoir in Portugal, from the Portuguese University of Évora, provided hourly data from 2017 and 2018. In addition, daily averaged values were also computed.

• In total, 27 Finnish lakes monitored by Finnish Environment Institute (SYKE) provided daily data (one measurement per day at 08:00 LT) from 2000 to 2016. In addition, summer month (June, July, and August) average values were calculated.

• Summer month average values were provided from the global inventory “Globally distributed LSWT collected in situ and by satellites; 1985–2009” (https://portal.edirepository.org/nis/mapbrowse?packageid=knb-lter-ntl.10001.3, last access: 29 August 2021). In total, there were 348 lakes over the globe with in situ data for 15 years (1995–2009).

Figure 5 provides the location of all the sources of in situ data used in this study. Considering the above observations as the truth, the MAE of ERA5 and ERA5-Land LSWT estimates was computed. The significance of these results was tested using the Kruskal–Wallis test by ranks. In addition, the bias distribution was also calculated and presented in 2-D bar graphs. It should be noted that the Alqueva reservoir and the 27 Finnish lake depths were verified (in situ vs. operational values) and are in good agreement; the lake depths provided by the “Globally distributed LSWT collected in situ and by satellites; 1985–2009” global inventory were only randomly verified (e.g. by comparison with scientific publications), which might add some uncertainty when interpreting the results.

Figure 5Location of lakes with in situ data used in this study; green dots are for hourly and daily data for the Alqueva reservoir (Portugal), red for daily and 3-summer-month averaged data (Finnish lakes), blue for the 3-summer-month average data for lakes all over the globe.

3.6 River discharge

The current version of CHTESSEL does not directly produce river discharge at the river basin scale. Instead, gridded surface and subsurface runoff from CHTESSEL is coupled to the LISFLOOD hydrological and channel routing model (Van Der Knijff et al., 2010). Coupling ERA5/ERA5-Land runoff with LISFLOOD allows for lateral connectivity of grid cells with runoff routed through the river channel to produce river discharge (m³ s⁻¹). This is the process used within the Global Flood Awareness System (GloFAS; https://www.globalfloods.eu/, last access: 29 August 2021). More details can be found in Harrigan et al. (2020). River discharge estimates from ERA5 (GloFAS-ERA5) and ERA5-Land (GloFAS-ERA5-Land) were obtained for the period January 2001 to December 2018. This is the common period for which reanalysis data were available at the time of this study. Estimates were resampled to the GloFAS 0.1^∘ gridded river network at a daily time step. As part of GloFAS, a database of global hydrological observations for 2042 stations is held, consisting predominantly (i.e. ∼ 75 %) of the Global Runoff Data Centre (GRDC) and supplemented by data collected through collaboration with GloFAS partners worldwide to improve spatial coverage. The locations of the stations have been matched to the corresponding cells on the 0.1^∘ GloFAS river network. Following Harrigan et al. (2020), a number of criteria were used to select stations for the evaluation:

• at least 4 years of data available between 2001 and 2018 (not necessarily contiguous),

• minimum upstream area of 500 km²,

• difference in catchment area supplied by the data provider and upstream area for the corresponding cell on the GloFAS river network must be within 20 %, and

• station with the longest record retained when multiple observation stations were matched to the same GloFAS river cell.

In addition to the above conditions, a first-order visual quality check on observed river discharge time series removed stations with erroneous data (for example, time series truncated above a threshold, showing several inhomogeneities or series monitoring an artificial canal instead of a river). This filtering procedure resulted in the selection of 1285 stations with drainage areas ranging between 575 and 4 664 200 km², and a median of 29 963 km². Following the methodology of Harrigan et al. (2020), hydrological performance was assessed using the modified Kling–Gupta efficiency (KGE^′) metric (Gupta et al., 2009; Kling et al., 2012). The KGE^′ is an overall summary measure consisting of three components important for assessing hydrological dynamics: temporal errors through correlation, bias errors, and variability errors:

where R is the Pearson correlation coefficient between reanalysis simulations (s) and observations (o), β is the bias ratio, γ is the variability ratio, μ the mean discharge, and σ the discharge standard deviation. The KGE^′ and its three decomposed components (correlation, bias ratio, and variability ratio) are all dimensionless, with an optimum value of 1. To evaluate the hydrological skill of GloFAS-ERA5-Land, the KGE^′ can be computed as a skill score, KGESS, with GloFAS-ERA5 used as the benchmark:

where KGEGloFAS−ERA5-Land′ is the KGE^′ value for the GloFAS-ERA5-Land reanalysis against observations, KGEGloFAS−ERA5′ is the KGE^′ value for the GloFAS-ERA5 benchmark against observations, and KGEperf′ is the value of KGE^′ for a perfect simulation, which is 1. KGESS=0 means the GloFAS-ERA5-Land reanalysis is no better than the GloFAS-ERA5 benchmark and thus has no skill, KGESS>0 indicates when GloFAS-ERA5-Land is considered skilful, and KGESS<0 is when the performance is worse than the GloFAS-ERA5 benchmark.

3.7 Energy fluxes

3.8 Skin temperature

The skin temperature is the theoretical temperature of the Earth's surface that is required to satisfy the surface energy balance. It represents the temperature of the uppermost surface layer, which has no heat capacity and thus can respond instantaneously to changes in surface fluxes. The top surface layer of the ECMWF's land surface model, CHTESSEL, covers the top 7 cm. In order to evaluate the skill of ERA5-Land land surface temperature (LST), NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) MYD11C3/MOD11C3 version 6 product was used in this study. It provides monthly LST and emissivity values for Aqua and Terra in a 0.05^∘ (5600 m at the Equator) latitude–longitude climate modelling grid (CMG), for day and night overpasses. Chen et al. (2017) recommended the use of the MODIS LST average ensemble (i.e. from Aqua-day, Aqua-night, Terra-day, and Terra-night) for climate studies and reported validation with 156 flux tower measurements (RMSE = 2.65, mean bias < ±1 K). In this study, the MODIS observational LST ensemble was constructed following Chen et al. (2017) and used as a reference for comparison with ERA-Interim, ERA5, and ERA5-Land monthly averaged LST data for the period January 2003 to December 2018 (16 years). MODIS data were firstly averaged to monthly timescales and then upscaled at two spatial resolutions: to 0.1^∘ for comparison to ERA5-Land, and to 0.25^∘ for comparison to ERA5 and ERA-Interim data. Bias, Pearson's correlation (R), and RMSE were computed on the full time series, and correlation was also computed on the anomalies (R_AN), calculated as departures from the monthly climatology.

4.1 Soil moisture

Figure 6 shows the evaluation results for ERA5-Land and ERA5 top soil moisture layer against 5 cm depth measurements for sites in Europe, Africa, and Australia. For sites in Europe, both ERA5-Land and ERA5 show similar STDD and bias distributions (Fig. 6a, c). The distribution of R is also similar but the median value obtained for ERA5-Land is slightly higher than for ERA5 (Fig. 6b). In contrast, ERA5 shows a slightly higher median for the anomaly correlation, although ERA5-Land shows a more compact distribution (Fig. 6d). In Africa, only 10 sensors at a depth of 5 cm were available. The STDD shows a larger distribution for ERA5-Land (Fig. 6e); however, the latter shows slightly lower bias and better R (Fig. 6f, g). The anomaly correlation is largely improved (Fig. 6h). Finally, in Australia, ERA5-Land box plots (Fig. 6i–l) clearly show lower STDD and bias, and higher R than ERA5 (both for the original time series and the anomaly times series). The evaluation obtained over Europe, Africa, and Australia shows an overall slightly better performance of ERA5-Land over ERA5; in particular, the anomaly correlation of the top layer is improved predominantly in the warmest climates. Nonetheless, these results for the top soil layer are not conclusive, partly because of the insufficient number of available stations.

Figure 6Box plots showing the evaluation of ERA5-Land and ERA5 top layer soil moisture against in situ measurements at 5 cm for sites in Europe (a–d), Africa (e–h), and Australia (i–l). Panels (a), (e), and (i) show the standard deviation of the difference (STDD), (b), (f), and (j) the Pearson correlation coefficient (R), (c), (g), and (k) the bias, and (d), (h), and (l) the Pearson correlation of the anomaly time series (R_AN). On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th (q₂₅) and 75th (q₇₅) percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers.

Download

To get more insight into the soil moisture evaluation, ERA5-Land and ERA5 were also evaluated over North America, where most of the in situ sensors are available, including several hundreds of sensors at 20 and 50 cm depths. In addition, the same evaluation was performed for ERA-Interim to address the evolution of ERA5 and ERA5-Land with respect to the previous reanalysis generation. Figure 7 shows the evaluation results for ERA5-Land, ERA5, and ERA-Interim soil moisture top three layers against measurements for sites in North America at 5 cm (Fig. 7a–d), 20 cm (Fig. 7e–h), and 50 cm depth (Fig. 7i–l). ERA5 and ERA5-Land surface soil moisture obtains quite similar results for all the metrics, although the median of the correlation and anomaly correlation are slightly better in ERA5-Land. ERA-Interim obtains lower values for these last two metrics, meaning that the land improvements in CHTESSEL compared to the old TESSEL scheme lead to better representation of the soil moisture dynamics. In deeper layers, the performance of ERA5-Land is better than that of ERA5, in particular for the third layer, for which ERA5-Land shows lower STDD (Fig. 7i) and higher correlation (Fig. 7j, l) than ERA5. It is worth mentioning that for the third soil layer of the model, ERA5 performs quite similarly to ERA-Interim, which is likely due to the initialization of the ERA5 streams with ERA-Interim soil moisture conditions.

Finally, the results were also analysed site per site taking into account the confidence interval of the Pearson correlation obtained for ERA5 and ERA5-Land versus the in situ measurements. A correlation difference is considered significant if the confidence intervals do not overlap. The Pearson correlation difference for ERA5 and ERA5-Land at 5 cm is significant for 382 sensors, of which 64 % shows higher correlation for ERA5-Land. The correlation difference at 20 cm is significant for 417 sites, of which 72 % show a higher correlation for ERA5-Land. Finally, the correlation difference at 50 cm is significant for 479 sites, of which 88 % show a higher correlation for ERA5-Land.

Figure 7Box plots showing the evaluation of ERA5-Land, ERA5, and ERA-Interim against in situ measurements at 5 (a–d), 20 (e–h), and 50 cm (i–l) over North America. Panels (a), (e), and (i) show the standard deviation of the difference (STDD), (b), (f), and (j) the Pearson correlation, (c), (g), and (k) the bias, and (d), (h), and (l) the Pearson correlation of the anomaly time series (R_AN). On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th (q₂₅) and 75th (q₇₅) percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers.

Download

4.2 Snow

Figure 8 shows the mean bias and RMSE of the SWE normalized by the mean value and standard deviation of the observations, respectively, for the 10 sites of Table 3. Each column displays the statistics of the nearest neighbour point, whereas the vertical bars represent the minimum and maximum values of the metrics computed over the four nearest points. Hence, this enables the characterization of the spatial variability of the errors given that many sites are located in complex terrains or coastal area.

Figure 8Summary statistics of the normalized mean bias (a) and RMSE (b) of snow water equivalent (SWE) at each ESM-SnowMIP site for ERA5-Land (red), ERA5 (green), and ERA-Interim (cyan). ESM-SnowMIP sites are described in Table 3. Each reanalysis box has a vertical line representing the variability at the four nearest grid points to the site location.

Download

Considering the mountain sites, ERA5-Land shows lower RMSE over the sites with moderate altitude, i.e. between 1300 and 2500 m (cdp, rme, wfj). ERA5-Land also clearly presents the lowest errors in the Arctic site (sod) and in the boreal forest sites (oas, obs, and ojp). Overall, the biases are smaller in ERA5-Land than in ERA5 (and ERA-Interim) at these sites. For the mountain sites, all reanalyses are characterized by a negative bias, which is very likely due to the smoothing of the orography at the resolution of the reanalysis. The higher horizontal resolution of ERA5-Land, compared to ERA5, helps to reduce the bias at cdp, rme, and wfj caused by a better orographic representation. However, the agreement with in situ observations at the sites located in very high mountains (snb and swa, located at altitudes greater than 3300 m) is slightly better with ERA5 than with ERA5-Land. It should be noted that the four nearest grid points indicate a much larger spread of errors for ERA5 than for ERA5-Land. Therefore, compensating errors could lead to better performance of ERA5 at the nearest grid point. The maritime site (sap) also shows lower RMSE in ERA5. At this site, the data assimilation can help in adding/removing snow mass for the right reason (snow density is overall well represented at this site; see time series in Fig. S1 in the Supplement), even though the spatial variability is higher in ERA5. On the contrary, at the forest sites (oas, obs, ojp), snow depth assimilation can remove snow mass to compensate for errors in snow density, the latter which are not considered in the assimilation system. Finally, noteworthy are the improvements in the transition between ERA-Interim and ERA5, in particular at Sodankylä (sod), which is caused by improved parameterizations of the snow model introduced between ERA-Interim and ERA5 productions (Dutra et al., 2010).

Figure 9 shows the maps of the RMSE difference between ERA5-Land and ERA5 snow depth estimates when compared to in situ observations of the GHCN network, for both North America and Europe. Over the US, and particularly over the Rockies region, ERA5-Land generally outperforms ERA5 in terms of lower RMSE (see Fig. 9a). In these highly complex terrain regions, the higher horizontal resolution of ERA5-Land adds value by providing more realistic orographic contours. However, over Europe (i.e. mainly Scandinavia; see Fig. 9b), where ERA5 uses a dense SYNOP network of observations in the snow assimilation system, ERA5 performs better than ERA5-Land.

Figure 9Snow depth RMSE difference between ERA5-Land and ERA5 (with respect to GHCN observations) for the months of December to June between 2010 and 2018. Negative values indicate better matching of ERA5-Land averaged snow depth with in situ measurements of the GHCN network; positive values indicate that ERA5 averaged snow depth estimates are matching better in situ measurements.

To further quantify the impact of the higher horizontal resolution of ERA5-Land on snow depth simulation, Fig. 10 shows the RMSE as a function of the height of each station. The stations were binned in height ranges of 250 m (0–250, 250–500, etc.) and the distribution of the RMSE is displayed for each height range using box plots. For heights below ≈ 1500 m a.s.l., ERA5 performs slightly better than ERA5-Land, while for heights between ≈ 1500 m a.s.l. and ≈ 3000 m a.s.l. ERA5-Land outperforms ERA5, as a result of the better resolution. For stations above 3300 m a.s.l., ERA5 performs better, which can be due to compensating errors, similarly to the ESM-SnowMIP sites. In addition, it should be noted that the number of sites at these very high altitudes is small, and therefore the statistical results should be interpreted with caution.

Figure 10Box plot of the snow depth RMSE distribution as function of the site altitude, for North America (30–80^∘ N, 60–160^∘ W, a) and Europe (25–80^∘ N, 10^∘ W–40^∘ E, b), for ERA5-Land (red), ERA5 (green), and ERA-Interim (cyan). Boxes extend between the lower (25 %) and upper quartiles (75 %), and the horizontal lines within each box represent the median value of the distribution. The black line represents the number of stations grouped at each bin.

Download

4.3 Lakes

The ERA5 and ERA5-Land hourly estimates of LSWT compared to the hourly data of the Alqueva reservoir showed that ERA5-Land obtained statistically significant MAE reduction by 2.2 % (see Table 4). While this is a positive result, ERA5-Land data showed some rapid changes in the lake mixed-layer depth (up to 5 m), which led to a quick rise of the lake surface water temperature (up to 20 ^∘C). The left plot of Fig. 11 shows the bias distribution for hourly data; even though the frequency of bias around zero is larger for ERA5-Land (indicating more accurate estimates of LSWT), errors due to unrealistic temperature rise are also visible at higher frequency between 3–4 ^∘C.

Table 4MAEs of the ERA5 and ERA5-Land LSWT estimates compared to in situ hourly, daily, and 3-summer-month average LSWT data. Hourly data correspond to data from the Alqueva reservoir in Portugal (2017 and 2018); daily data are from the Finnish lakes (2000–2016) and averaged hourly data from the Alqueva reservoir; summer month average data are from a global inventory (1985–2009) and the Finnish lakes. The first column is for the data type (hourly/daily/summer month average), second column provides the number of lakes for each category, third column shows the MAE (^∘C) of ERA5 LSWT estimates compared to in situ measurements, and fourth column is the same as third column but for ERA5-Land.

Download Print Version | Download XLSX

Figure 11LSWT bias (observations minus ERA5 or ERA5-Land estimates) distribution (^∘C) for hourly data from the Alqueva reservoir (a) and daily data from the Alqueva reservoir and the Finnish lakes (b). The dark red colour represents the area where both ERA5 (pink colour) and ERA5-Land (violet colour) bias distribution overlap.

Download

The comparison of the ERA5 and ERA5-Land daily averaged LSWT with Alqueva reservoir and the 27 Finnish lakes data also showed statistically significant reduction of ERA5-Land MAE (among 28 lakes) by 1.2 % (from 2.71 to 2.68 ^∘C; see Table 4). It should be noted that the lake parametrization scheme relies heavily on the lake depth input data. With the increase of resolution in ERA5-Land, lake depth values change as well. For only those lakes whose depth became more realistic in ERA5-Land as result of the higher resolution, the MAE was statistically significantly reduced by 23.8 %; however, for only those lakes whose depth remained unchanged despite the change of resolution, the MAE was still statistically significantly reduced by 9.6 %, showing the positive influence of higher-resolution atmospheric input. The right plot of Fig. 11 shows, for daily data, a reduced amount of large positive errors (cold bias) in ERA5-Land and a larger frequency of errors around zero bias. Finally, only 272 lakes with averaged in situ measurements for the summer months, with 10- to 15-year data records, could be used for comparison with ERA5 and ERA5-Land LSWT estimates. The comparison showed statistically significant MAE increase by 1.0 % for ERA5-Land compared to ERA5 (see Table 4). If one bears in mind that FLake was not designed for saline, fed by glaciers or mountains and warm lakes (that excludes 26 exceptional lakes) the MAE of the remaining 246 lakes is similar in ERA5 and ERA5-Land. The right plot of Fig. 12 shows the geographical distribution of the 26 exceptional lakes and their averaged ERA5-Land MAE, which shows large LSWT errors, even beyond 10 ^∘C. Note that errors for these lakes are very similar for ERA5. If, as a result of the higher resolution, only lakes whose depth became more realistic in ERA5-Land compared to ERA5 are taken into account, the MAE is statistically significantly reduced by 1.4 % in ERA5-Land. FLake was specially designed for medium depth (under 50 m) lakes. Several lakes from this database are however very deep. Nevertheless, computing the statistics for only medium depth lakes did not show any extra improvement in ERA5-Land data. The left plot of Fig. 12 presents the bias distribution of non-exceptional lakes (246) that shows wider error distribution for both ERA5 and ERA5-Land compared to hourly and daily distributions.

Figure 12LSWT bias (observations minus ERA5 or ERA5-Land estimates) distribution (^∘C) for the 3-summer-month average data based on non-exceptional lakes (a), and geographical distribution of the 26 exceptional lakes (glacier-fed, saline, and warm lakes) ERA5-Land MAE (b) with colour bar in ^∘C.

4.4 River discharge

Results of river discharge performance for GloFAS forced with ERA5 and ERA5-Land are shown in Fig. 13. The overall global median KGE^′ across 1285 observation stations improves from 0.26 (with an interquartile range of −0.04 to 0.49) for GloFAS-ERA5 to 0.37 (0.08, 0.57) for GloFAS-ERA5-Land (Fig. 13a). In the decomposition of the KGE^′, the global median Pearson correlation also improves from 0.60 (0.43, 0.75) to 0.64 (0.49, 0.77) (Fig. 13b), and global median bias ratio improves from 0.73 (0.50, 1.13) to 0.89 (0.66, 1.15) (Fig. 13c), which is equivalent to a 16 % reduction in overall bias. There is very little difference in variability errors between GloFAS-ERA5 and GloFAS-ERA5-Land (Fig. 13d).

Figure 13Cumulative distribution functions (CDFs) of hydrological performance metrics across 1285 observation stations. Modified Kling–Gupta efficiency (KGE^′) (a) with decomposition of KGE^′ into Pearson correlation (b), bias ratio (c), and variability ratio (d) for GloFAS-ERA5-Land (red line) and GloFAS-ERA5 (blue line). The red dot marks the optimum value for each metric.

Download

Figure 14 shows the spatial distribution of river discharge skill, measured by the KGESS. GloFAS-ERA5-Land shows positive skill compared to the GloFAS-ERA5 benchmark in 65 % of stations with a global median KGESS of 0.08 (−0.06, 0.25). Largest improvements in skill are found in North America, Europe, northern Russia, southern Africa, and Australian catchments. There is a substantial decrease in skill (i.e. KGESS < −0.2) when forcing GloFAS with ERA5-Land runoff instead of ERA5 in 11 % of stations, mainly located in the western US and South America. Care must be taken in spatial representativeness of these results, as the observation network is sparse in some regions of the world, particularly in large parts of Africa and Asia.

Figure 14Modified KGESS for GloFAS-ERA5-Land river discharge reanalysis against the GloFAS-ERA5 benchmark across 1285 observation stations. Optimum value of KGESS is 1. Blue (red) dots show catchments with positive (negative) skill.

4.5 Energy fluxes

4.5.1 Evaluation against eddy-covariance sites

The left panel of Fig. 15 shows the violin plots of the ERA-Interim and ERA5-Land turbulent fluxes (compared to in situ eddy-covariance measurements).

Figure 15Violin plots showing the temporally and spatially averaged statistics of the surface latent heat flux (λρE), surface sensible heat flux (H), and Bowen ratio (β) from (a) ERA-Interim (green) and ERA5-Land (yellow) and (b) ERA5 (green) and ERA5-Land (yellow). Statistics are calculated with respect to in situ eddy-covariance measurements at both 3-hourly and daily (24 h) temporal resolutions. Violin plots represent the distribution of the individual validation statistics with indication of the median and interquartile range, and are calculated using a kernel density estimation approach. Statistics include the bias, standardized MAE, and anomaly correlation coefficient (R_AN). The scale used to show the bias distribution of β is that of the right y axis at the top row of each panel.

Download

ERA5-Land compares systematically better to in situ measurements than ERA-Interim, with higher R_AN values for the heat fluxes and the Bowen ratio and for both subdaily and daily temporal resolutions, as well as presenting lower bias and MAE. This result is expected, as the land surface model used in ERA5-Land benefits from many improvements compared to that of ERA-Interim (see Sect. 2.4). The right panel of Fig. 15 shows the same plot but comparing ERA5 and ERA5-Land. The violin plots are much more similar, biases for both fluxes are only marginally better in ERA5-Land (the median of the distribution is nearly the same, but the 75 % percentile is always lower), and MAE is typically lower for λρE in ERA5-Land (except at 1-hourly resolution) but higher for H. On average, correlations for ERA5-Land are only better for λρE at daily resolution and for the Bowen ratio. It also should be noted that all statistics are on average better for H than for λρE, both in ERA5 and in ERA5-Land (with lower bias, lower MAE, and higher R_AN); this matches the results reported by Balsamo et al. (2015) and Martens et al. (2020). To investigate potential areas of significant improvement/degradation, the left panel of Fig. 16 shows the maps of R_AN differences of λρE between ERA5-Land and ERA5, over the continental United States (CONUS). R_AN is typically better for ERA5-Land across most US stations, especially near the coasts and around the lakes, where high resolution is more important. In Europe (Fig. 16, middle panel), results are very mixed, but for most of the stations in the Alps ERA5-Land performs worse than ERA5. Although for Australia ERA5-Land is slightly better than ERA5 (see right panel of Fig. 16), it is not the case for all stations. The results for the Bowen ratio are aligned with those of λρE (see Fig. S3), whereas those of H are more favourable for ERA5 (see Fig. S2).

Figure 16Standardized anomaly correlation (R_AN) difference of the surface latent heat flux between ERA5-Land and ERA5 with respect to eddy-covariance data. Blue colours indicate that the anomaly correlation of ERA5-Land with respect to eddy-covariance measurements is higher than for ERA5, whereas red colours indicate the opposite.

To study the influence of the eddy-covariance sites' altitude and air temperature in the heat fluxes, Fig. 17 shows the standardized mean R_AN (circles) and MAE (squares) of H and λρE as a function of the station altitudes and air temperature. The results are in line with those presented in the violin plots. ERA5-Land clearly performs better for λρE and the Bowen ratio, and the largest differences between ERA5 and ERA5-Land are obtained for the sites located at high altitudes, where more extreme values are obtained and forcing errors are increased. Although for H the differences are smaller, ERA5 performs overall better than ERA5-Land and the strongest differences are also seen at high-altitude sites.

Figure 17Standardized anomaly correlation difference (circles) and standardized mean absolute differences (squares) between ERA5-Land and ERA5 latent heat flux (a), sensible heat flux (b), and Bowen ratio (c), grouped as a function of stations temperature and altitude. The size of circles and squares is proportional to the number of eddy-covariance towers. Green values denote better matching of ERA5-Land with in situ data.

Download

4.5.2 Evaluation using GLEAM

The turbulent fluxes estimated from GLEAM forced with temperature and surface net radiation from ERA5 and ERA5-Land were also computed. ERA5-Land surface fluxes perform better than GLEAM + ERA5-Land when confronted to in situ measurements (see the violin plots in Fig. S4), except in terms of bias, in agreement with Martens et al. (2020). This result reflects the quality of the CHTESSEL land surface model, which has greatly improved the parametrization of turbulent energy fluxes and evidences the added value compared to a simpler model as in GLEAM, designed to consider only input variables that can be observed from satellite. Figure 18 is similar to Fig. 17 but comparing ERA5-Land and GLEAM + ERA5-Land, as a function of elevation and air temperature, at daily resolution. As expected, almost for all stations temperature and elevation, ERA5-Land heat fluxes match better the in situ reference measurements. The exceptions are towers at the highest altitudes where forcing errors are generally larger. Comparing GLEAM + ERA5 and GLEAM + ERA5-Land, the differences are very small, suggesting that the near-surface air temperature and net radiation are quite similar in both reanalyses for the locations of the eddy-covariance sites.

Figure 18Standardized anomaly correlation difference (circles) and standardized mean absolute differences (squares) between ERA5-Land and GLEAM + ERA5-Land latent heat flux (a), sensible heat flux (b), and Bowen ratio (c), grouped as a function of stations' temperature and altitude. The size of circles and squares is proportional to the number of eddy-covariance towers. Green values denote better matching of ERA5-Land with in situ data.

Download

4.6 Skin temperature

Figure 19 shows global maps of ERA5-Land and ERA-Interim mean LST for the time period 2003–2018, as well as global maps of their correlation and RMSE with respect to the MODIS LST average ensemble, constructed as indicated in Sect. 3.8. While differences in mean LST are apparently minimal among the two simulated products, there is a better correlation of ERA5-Land with MODIS in the tropical band. The low correlation obtained between the simulated and the remotely sensed LST over tropical forests should be taken with caution since persistent cloud cover as well as cloud contamination effects are known to have an impact on remotely sensed observations, leading to important differences in their absolute values (Jiménez-Muñoz et al., 2016; Gomis-Cebolla et al., 2018). Notably, ERA5-Land systematically compares better than ERA-Interim to MODIS LST in terms of RMSE, particularly at high latitudes. As shown in the averaged statistical scores on Table 5, there is an overall improvement of ERA5 and ERA5-Land LST simulations with respect to ERA-Interim, yet only modest improvements are achieved in ERA5-Land compared to ERA5 LST, mainly on bias and RMSE. The better correlation of MODIS with ERA5 and ERA5-Land than with ERA-Interim is more evident in R_AN maps (see Table 5 and Fig. S5).

Figure 19Global statistics of ERA5-Land (a, c, e) and ERA-Interim (b, d, f) LST (in K) for the time period 2003–2018. Panels (a) and (b) are the global maps of mean LST, (c) and (d) are the correlation maps, and (e) and (f) are the RMSE maps, with respect to the MODIS LST average ensemble. Note that the LST here refers to the skin temperature.

To further scrutinize the impact of the ERA5-Land higher horizontal resolution on LST simulation, ERA5 LST was resampled to 0.1^∘ and then global maps of the RMSE difference (ΔRMSE) between ERA5 and ERA5-Land with respect to the MODIS LST average ensemble were calculated. The improvement of ERA5-Land on LST in coastal areas is clearly visible in the ΔRMSE maps, as shown in Fig. 20 for the European domain. The time series at the coastal pixel in Norway illustrates how ERA5-Land LST simulations compare better to the observational record by capturing more closely the annual cycle in its whole range of variability, and therefore decreasing bias and RMSE, while the correlation with the MODIS product is very good and almost similar for both reanalyses. The time series at the other two specific locations illustrate regions of complex topography exhibiting slightly lower RMSE of ERA5-Land (Iceland), and of ERA5 (Alps), with respect to the MODIS LST average ensemble. These results highlight improved LST simulations of ERA5-Land on coastal regions.

Figure 20Map of RMSE differences (ΔRMSE (K)) between ERA5 and ERA5-Land with respect to the MODIS LST average ensemble for the time period 2003–2018. Red values denote lower RMSE of ERA5-Land with respect to the observational record. Time series of ERA5, ERA5-Land, and MODIS time series at selected pixels: coast of Norway (61.6^∘ N, 5.3^∘ E), the Alps (47.4^∘ N, 11.15^∘ E), and Iceland (64.5^∘ N, 20.26^∘ W). The differences in RMSE (ΔRMSE) and R (ΔR) between ERA5 and ERA5-Land with respect to MODIS LST at each location are reported.

Table 5Spatially and temporally averaged bias, RMSE, R, and anomaly correlation coefficient (R_AN) between ERA-Interim, ERA5, and ERA5-Land estimates and the MODIS land surface temperature average ensemble for the time period 2003–2018.

Download Print Version | Download XLSX