ERA5-Land: A state-of-the-art global reanalysis dataset for land applications

Framed within the Copernicus Climate Change Service (C3S) of the European Commission, the European Centre for Medium-Range Weather Forecasts (ECMWF) is producing an enhanced global dataset for the land component of the 5 generation of European ReAnalysis (ERA5), hereafter named as ERA5-Land. Once completed, the period covered will span from 1950 to present, with continuous updates to support land monitoring applications. ERA5-Land describes the evolution of the water and energy cycles over land in a consistent manner over the production period, which among others could be used 5 to analyse trends and anomalies. This is achieved through global high resolution numerical integrations of the ECMWF land surface model driven by the downscaled meteorological forcing from the ERA5 climate reanalysis, including an elevation correction for the thermodynamic near-surface state. ERA5-Land shares with ERA5 most of the parametrizations that guarantees the use of the state-of-the-art land surface modeling applied to Numerical Weather Prediction (NWP) models. A main advantage of ERA5-Land compared to ERA5 and the older ERA-Interim is the horizontal resolution, which is enhanced globally 10 to 9 km compared to 31 km (ERA5) or 80 km (ERA-Interim), whereas the temporal resolution is hourly as in ERA5. Evaluation against independent in situ observations and global model or satellite-based reference datasets shows the added value of ERA5-Land in the description of the hydrological cycle, in particular with enhanced soil moisture and lake description, and an overall better agreement of river discharge estimations with available observations. However, ERA5-Land snow depth fields present a mixed performance when compared to those of ERA5, depending on geographical location and altitude. The 15 description of the energy cycle shows comparable results with ERA5. Nevertheless, ERA5-Land reduces the global averaged root mean square error of the skin temperature, taking as reference MODIS data, mainly due to the contribution of coastal points where spatial resolution is important. Since January 2020, the ERA5-Land period available extends from January 1981 to near present, with 2 to 3 months delay with respect to real-time. The segment prior to 1981 is in production, aiming to a

ERA5-Land produces a total of 50 variables describing the water and energy cycles over land, globally, hourly and at a spatial 85 resolution of 9 km, matching the ECMWF triangular-cubic-octahedral (TCo1279) operational grid (Malardel et al., 2016). For a full list of the available fields in the ERA5-Land catalogue, see Table A2 in the appendix. The production is conducted in three segments or streams. The reason is twofold: 1) it allows the production of parallel streams, therefore accelerating the production and the public availability of the data, 2) the atmospheric forcing necessary to produce ERA5-Land is derived from ERA5, thus the production needs the corresponding segment of ERA5 completed for the same time period. Figure  year spin-up (1949). The red diamond presents the end of each stream.
assimilate observations directly. The observations influence the land surface evolution via the atmospheric forcing. Forcing air temperature, humidity and pressure are corrected using a daily lapse rate derived from ERA5. After that, the land surface 95 model is integrated in 24 h cycles providing the evolution of the land surface state and associated water and energy fluxes. In addition to the hourly data, monthly means are also computed (Muñoz-Sabater, J., 2019b). Figure 2 shows a diagram of the algorithm used for each 24 h production cycle. The most important components of the production algorithm are presented in the following subsections.

100
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 Figure 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, Figure 3 shows time series of deep soil moisture in a band of latitude between 60°S 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 105 ERA-Interim soil moisture initial conditions, which has a different climatology than ERA5. In ERA5 one-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, what leads to discontinuities between two production streams, as shown in Fig. 3. 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 three further spin-up years to allow a long spin-up period (see Fig. 1). While this strategy 110 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 st of 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 115 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 the 1 st of January for the period 2001-2018, and then allowing the system to spin-up for three 120 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 one spin-up year was feasible. The latter was limited by the availability of forcing data.

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 type, and vegetation type and vegetation cover. In addition surface albedo and leaf area index 125 are prescribed as monthly climatologies. The complete list of time invariant fields with their information source is provided in table A1 of the appendix.

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: 130 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 all the production period, and is the result of the assimilation of a large number of conventional meteorological and satellite observations through a 4-Dimensional Variational assimilation 135 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 to (2) reduce dependencies on external data that would limit the near-real time data availability. However, air temperature, humidity and 140 pressure are corrected for the altitude differences between ERA5 and ERA5-land grids. This correction involves 4 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.  Primary production (GPP) and the Ecosystem respiration (Reco)) in a modular way with the Jarvis approach that computes the 165 stomatal conductance without affecting the transpiration components (Boussetta et al., 2013a). 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 175 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.

Data and Evaluation strategy
To evaluate the quality of the ERA5-Land fields, several key variables of the water and energy cycles were selected and compared to available in situ observations and to a series of reference datasets. Note that the list of evaluated variables and reference datasets is not exhaustive and was based on factors such as availability of data at the time of the evaluation. ERA-Interim and ERA5 reanalyses were included in the comparison, aiming at showcasing the progress of operational reanalyses 185 at ECMWF. The variables evaluated are soil moisture, snow depth, lake surface water temperature and river discharge for the water cycle, the sensible and latent heat fluxes (the latter also a component of the water cycle), the Bowen ratio and skin temperature for the energy cycle. This section describes the supporting datasets used in the evaluation exercise and the metrics employed to assess the quality of ERA5-Land fields.

190
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 still today is widely used by the scientific community. It includes information on multiple land and atmospheric variables, and it is available from January 1979 to 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 hours and 3 hours 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 as 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 200 assimilates SYNOP 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).

ERA5
ERA5 is the latest comprehensive ECMWF reanalysis and has replaced ERA-Interim. It is based on a version of the ECMWF 205 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 days. 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 utilise the output from the 210 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 CMIP-5 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 land data assimilation system (LDAS). The

215
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 1m of soil, and a one-dimensional Optimal Interpolation for soil, ice and snow temperature.
Details of the ERA5 configuration are described in Hersbach et al. (2020) that 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 220 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%3A+data+documentation).
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 Sore (CDS). ERA5 and ERA5-Land have 225 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.

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 230 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(Dorigo et al., , 2021. Three reanalysis soil layers 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 was also evaluated. In addition, ERA5, ERA5-Land and ERA-Interim soil moisture for the second (7-28 cm) and third layer (28-100 cm) were evaluated against in situ measurements at 20 cm 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 240 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 days 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, and equivalent to the unbiased RMSD), the bias and the Pearson correlation coefficient (R). The latter was also computed for the anomalies 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 q 75 + 1.5 × (q 75 − q 25 ) or less than q 25 − 1.5 × (q 75 − q 25 ), with q 25 and q 75 the 25 th and 75 th percentiles, respectively.

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 255 cover and on different types of climates. Table 3   reported is lower than 1 cm in more than 5% of the total number of days were removed. Finally, stations located on coastal areas with more than 50% of water in the pixel and on permanent snow area (glaciers) were removed. More than 6000 stations passed the quality filters, and their location can be consulted 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.

Lakes
In 2015, the CHTESSEL land surface scheme of the operational IFS introduced a lake tile, which represents lakes, reservoirs, rivers and coastal (sub-grid) waters, and is based on the FLake (Fresh-water Lake) model of Mironov et al. (2010). 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 275 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. 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 280 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: 13 -The Alqueva reservoir in Portugal, from the Portuguese University of Evora, providing hourly data from 2017 and 2018.
In addition, daily averaged values were also computed,  for daily and three summer month averaged data (Finnish lakes), blue for three summer month average data for lakes all over the globe.

River discharge
The current version of CHTESSEL does not directly produce river discharge at the river basin scale. Instead, gridded surface and sub-surface runoff from CHTESSEL is coupled to the LISFLOOD hydrological and channel routing model (Van Der Kni-300 jff 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 3 · s −1 ). This is the process used within the Global Flood Awareness System (GloFAS; https://www.globalfloods.eu/). 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 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: 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: is the KGE value for the GloFAS-ERA5 benchmark against observations, and KGE perf is the value of KGE for a perfect simulation which is 1. A KGESS = 0 means the GloFAS-ERA5-Land reanalysis is no better than the GloFAS-ERA5 benchmark so has no skill, KGESS > 0 when GloFAS-ERA5-Land is considered skilful, and KGESS < 0 when the performance 335 is worse than the GloFAS-ERA5 benchmark.

FLUXNET data
The evaluation of the ERA5-Land turbulent fluxes estimates was conducted mostly following the method of Martens et al.
(2020). Surface sensible and latent heat fluxes (also denoted in this paper as H and λρE, respectively) derived from ERA 340 reanalyses, as well as their ratio (i.e. the Bowen ratio, hereafter denoted as β), were compared to measurements from the FLUXNET 2015 synthesis data set (Pastorello et al., 2020). The period under evaluation was based on the availability of reanalysis data at the time of the comparison, and therefore, unlike in Martens et al. (2020), the evaluation period is constrained to 2001-2014. Following Martens et al. (2017), the in situ flux data was subjected to quality control, including (1) the removal of rainy intervals, during which eddy-covariance measurements are typically unreliable, and (2) the removal of gap-filled 345 records to retain only the actual measurements from the eddy-covariance sites. After quality control only sites with a minimum record of five years were retained too. In total, 65 eddy-covariance sites remained after quality control and were used as in situ reference data. Note that the measured energy fluxes used as reference in this paper were not corrected for energy balance closure because the number of towers used for validation would be drastically reduced, as the ground-heat flux is also needed and is not available from many towers. Some authors have already highlighted the lack of closure in the energy balance at 350 eddy-covariance sites and a consequential tendency to underestimate the latent heat flux (Wilson et al., 2002;Ershadi et al., 2014;Jiménez et al., 2018). The reference sites are mainly distributed across the continental US, Europe and Australia (see their location in Fig. 17a, b and c, respectively). For each eddy-covariance site, the in situ measurements were aggregated from their native temporal resolution to hourly, 3-hourly and daily intervals. In addition, standardized anomalies were calculated by subtracting for each time interval the climatological expectation (i.e. the average value across the entire record for that interval) 355 and dividing by the standard deviation of that climatology. While the comparison to raw time series may mask the influence of short-term meteorological anomalies on surface energy partitioning (as the temporal variability of turbulent fluxes typically depends strongly on the seasonality of its main drivers), the comparison to anomaly time series reflects the response to short-term meteorological conditions. The Bowen ratio was only calculated at daily temporal resolution for numerical instability reasons. As described in Martens et al. (2020), outliers in the time series of the Bowen ratio, for both the reanalyses and in situ 360 data, were masked using a quantile-based approach. The bias, i.e., the difference between raw in situ time series and reanalysis estimates, the standardized MAE and the anomaly Pearson correlation coefficient (R AN ) were computed. The results are shown in the form of violin plots in section 4.5.1.

GLEAM
The Global Land Evaporation Amsterdam Model (GLEAM; Miralles et al. (2011), Martens et al. (2017) is used in this study 365 with two objectives: a) to compare the GLEAM evaporation estimates directly to the ERA5 and ERA5-Land estimates, which in turn will assess the skill of the underlying land-surface model to simulate turbulent heat fluxes, and b) as an intermediate tool to assess the quality differences of key input meteorological drivers of the turbulent fluxes computed by GLEAM. GLEAM is a process-based, yet semi-empirical, model that computes total evaporation and its separate components over continental masses at global scale. A detailed description of this model can be found in Martens et al. (2017) and Miralles et al. (2010, 370 2011). In this paper the version 3 (v3) of the GLEAM algorithm is used, and forced with the same database as in the official v3.4a dataset (see www.gleam.eu), including near surface air temperature and surface net radiation from ERA5; this dataset is hereafter referred to as GLEAM+ERA5. Likewise, a version of GLEAM run with ERA5-Land air temperature and surface net radiation is referred to as GLEAM+ERA5-Land. Analogous comparisons between GLEAM+ERA5 and GLEAM+ERA-Interim (using forcing fields from ERA-Interim) can be found in Martens et al. (2020). Surface latent heat flux, surface sensible 375 heat flux and the Bowen ratio from GLEAM+ERA5 and GLEAM+ERA5-Land were evaluated against the eddy-covariance data described in section 3.7.1 at daily time scales.

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 so can respond instantaneously  (Figs. 6a, c). The distribution of the 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.   400 6f,g). The anomaly correlation is largely improved (Fig. 6h). Finally, in Australia ERA5-Land boxplots (Figs. 6i-l) clearly show lower STDD and bias, and higher R than ERA5 (both for the original time series and the anomalies 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.

405
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     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, enabling the characterisation of the spatial variability of the errors given that many sites are located in complex terrains or coastal area.
Considering the mountain sites, ERA5-Land shows lower RMSE over the sites with moderate altitude, i.e, between 1300 m Snow water equivalent era5land era5 erai 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 4 nearest grid points to the site location and 2500 m (cdp, rme, wfj). ERA5-Land also clearly presents the lowest errors in the Arctic site (sod) and in the Boreal forest 430 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 reducing 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 435 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 performances 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 supplementary Fig. S1), even though the spatial variability is higher in ERA5. On the contrary, at the forest sites (oas, obs, ojp), snow depth assimilation can 440 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 Sodankyla (sod), which is caused by improved parametrizations of the snow model introduced between ERA-Interim and ERA5 productions . Figure 9 shows the maps of the RMSE difference between ERA5-Land and ERA5 snow depth estimates when compared 445 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. 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 boxplots. 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 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.

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-

460
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). 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 errors.
The comparison of the ERA5 and ERA5-Land daily averaged LSWT with Alqueva reservoir and the 27 Finnish lakes data 465 also showed statistically significant reduction of ERA5-Land MAE (among 28 lakes) by 1.2 % (from 2.71 • C 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 470 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-15 years-long 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 475 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 480 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 errors distribution for both ERA5 and ERA5-Land compared to hourly and daily distributions.
25 Figure 11. Lake Surface Water Temperature bias (observations minus ERA5 or ERA5-Land estimates) distribution ( • C) for hourly data from the Alqueva reservoir (left) and daily data from the Alqueva and the Finnish lakes (right). The dark red colour represents the area where both ERA5 (pink colour) and ERA5-Land (violet colour) bias distributions overlap.  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 495 ERA5 in 11 % of stations, mainly located in 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. 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), 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 To study the influence of the eddy-covariance sites altitude and air temperature in the heat fluxes, Fig 17 shows  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.

Evaluation using GLEAM
The turbulent fluxes estimated from GLEAM forced with temperature and surface net radiation from ERA5 and ERA5-Land   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 545 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 in the supplementary material).
To further scrutinize the impact of the ERA5-Land higher horizontal resolution on LST simulation, ERA5 LST was resam-

Discussion and conclusions
This paper presents the new global ERA5-Land reanalysis. When the historical part is completed it will provide a detailed record of the land surface evolution from 1950 to present through a series of key land surface variables representing the water and energy cycles. The temporal resolution is hourly and the horizontal resolution is 9 km, making it unique and suitable for   scales. ERA5 root-zone soil moisture is penalized by initialization from ERA-Interim, which leads to a multi-annual artificial trend of the time series over arid to semi-arid climates, where the standard deviation of soil moisture is low. On the contrary, the shallow surface layer responds quickly to short-term meteorological forcing variables, and hence, for the top layer, ERA5-Land only improves slightly on ERA5 estimates at the reference stations. The main added value of high resolution at the top layer is to account for the correct soil type that changes the saturation level and the minimum level at which evapotranspiration 575 no longer occurs. Consequently, ERA-Interim, with coarser resolution, shows the lowest performance of the three reanalyses datasets. Supporting these findings, ERA5-Land participated in a soil moisture intercomparison of 18 products against large number (826) of in situ stations and performed strongly (Beck et al., 2021).
(b) Snow: ERA5-Land snow mass and snow depth estimates are improved on mid-altitude mountains, where the higher res- olution provides better orographic contours and the temperature correction is also important. However, ERA5 snow depth 580 estimates match measurements slightly better at the highest mountains (>3300 m). At these altitudes, uncertainties in both the forcing and in the model parametrization are larger and can contribute to the growth of errors; processes like snow transport and sublimation, which are not considered in the current snow scheme formulation, can be important at these altitudes, as well as errors in the amount of solid precipitation. It is also worth noting that for high altitude stations, the spread of errors is larger in ERA5 compared to ERA5-Land, which could lead to error compensation. The results per continent show the best performance 585 of ERA5-Land snow fields in the US, especially in complex orographic areas, where spatial resolution is foremost. It implies a more accurate near surface air temperature, especially over mountains, which leads to a better representation of the sensible heat flux and therefore better estimates of the snow depth. Also, there is a limited number of snow depth in situ observations assimilated in ERA5 over the US (de Rosnay et al., 2015), and therefore the assimilation system hardly compensates for systematic model/snowfall biases. Contrarily, ERA5 performs better over Scandinavian countries. In general snowfall in the IFS 590 suffers from a systematic positive bias (snow overestimation) in Europe. Although data assimilation is intended to correct for random errors, snow depth reports assimilated in ERA5 alleviate the systematic snow excess where these reports are assimilated. The improvements of the snow scheme over the previous decade is also reflected in the overall better scores of ERA5 compared to ERA-Interim.
(c) Lakes: LSWT is slightly improved in ERA5-Land at hourly and daily averaged steps with respect to in situ data at sites in 595 Portugal and Finland. Note however that at punctual times of the warm period in Alqueva, sudden jumps of the ERA5-Land mixed layer depth leads to unrealistic LSWT. For instance, on 3 August 2018, from 8 am to 9 am local time, the mixed layer depth jumped 5.3 m, leading to a 17.6°C increase and increasing the Alqueva reservoir LSWT to 40.7°C. Jumps are caused by a simplistic parametrization of the summer stratification in the lake scheme, as temperature is computed differently in calm and turbulent conditions. Ways to limit these jumps and to enable a smoother and more realistic transition are currently being 600 investigated at ECMWF. The lake parametrization scheme relies heavily on accurate lake depth input data, and therefore poor data impacts its ability to accurately simulate lake parameters such as LSWT. The higher resolution of ERA5-Land allows, in many cases, a more accurate specification of lake depth. Thus, isolating lakes whose depth is more realistic using the ERA5-Land grid (i.e. only lakes whose depth in the ERA5-Land grid matches better in situ observations), the LSWT MAE is reduced by more than 20%. The positive influence of a higher resolution atmospheric forcing was also verified by isolating lakes whose 605 depth remained unchanged either using ERA5 or ERA5-Land.
Both reanalysis LSWT estimates were also compared using as reference a global inventory (1995-2009) based on summer observations from satellite sensors. In this case, the performance of both reanalysis is quite similar (excluding exceptional lakes where Flake performs poorly). The uncertainty of lake depths in the global inventory is larger and it could have an impact on these results. However, even when comparing only lakes with a depth less than 50 m (Flake was designed for medium-depth 610 lakes) the performance of ERA5-Land LSWT does not show significant improvement compared to ERA5 (based on MAE). This might be due to the averaging technique used for the in situ measurements, which uses satellite measurements representing one instant in time rather than continuous hourly data. Nevertheless, in summer months, ERA5 LSWT biases are on average 2.2°C cooler than observations, whereas ERA5-Land LSWT is just 1.3°C cooler.
(d) River discharge: River discharge integrates the various components of the water cycle. With improved soil moisture, snow 615 and lakes characterization in ERA5-Land, all these being physical inputs of river discharge generation, it is not surprising that the GloFAS river discharge obtained overall better scores when used with input data from ERA5-Land (GloFAS-ERA5-Land).
This improvement in the skill is likely due to a number factors, such as higher horizontal resolution (from which smaller catchments benefit), improved representation of root-zone soil moisture and temperature (which contributes to increase the correlation and reduce biases with river discharge observations) and snow processes (which should be an advantage in northern 620 latitude catchments). Also, the fact that ERA5-Land does not directly assimilate any observation, as is the case for ERA5, has shown to have positive impacts on the closure of the water balance (Zsoter et al., 2019). They have shown that the snow data assimilation is detrimental to the hydrology in large parts of the snow impacted Northern Hemisphere in the ERA5 experiment, compared with the offline simulation without coupling and land data assimilation (such as ERA5-Land). The average snow increments in ERA5 are negative, due mainly to the too slow snow melt in CHTESSEL, which consequently removes water from 625 the hydrological system. This contributes to decreasing the dominantly negative biases in a large area with subsequently deteriorating the hydrological performance. The exceptions are in the West US and Amazonian basins. For the former, a later and larger snow melt due to an excess of accumulated snow over high peaks could lead to a decrease in correlation and increased biases with respect to in situ observations. The degradation in the Amazonian basin is small, linked to an underestimation of the surface and subsurface runoff (not shown) which in turn is likely caused by a slight overestimation of evaporation. (f) Skin Temperature: The skin temperature (in this paper referred as LST), as a variable reacting quickly to any change in surface fluxes, shows modest improvements in ERA5-Land compared to ERA5. However, as it has been the case for all variables 650 evaluated in this paper, ERA5 and ERA5-Land LST are improved when compared to ERA-Interim. ERA5-Land obtains the lowest global averaged RMSE with respect to MODIS data, which is partly due to the contribution of the coastal points that simulate better the amplitude of the annual cycle of LST, and is a consequence of the higher spatial resolution of ERA5-Land.
Other small disagreements between ERA5 and ERA5-Land are found over complex terrains, but they do not seem to favour any particular reanalysis.
With the above results, one can conclude that the horizontal resolution matters and is a very important aspect in the accurate simulation of the spatial and temporal evolution of the hydrological cycle. However, this paper could only provide evidence of a modest improvement of the surface fluxes of ERA5-Land compared to ERA5. The latter conclusion is based on an evaluation with respect to a small number of available samples. Other important aspects of the added value of ERA5-Land 660 are the production speed (that allows cutting-edge land surface modeling advances to be incorporated more rapidly) and the consistency presented over multi-decadal time scales (that could set the basis to enable reliable trend analyses), all of them making ERA5-Land a state-of-the-art dataset for multiple land applications. The reduced impact of discontinuities by using longer spin-up periods, could also be a crucial factor to obtain accurate trends over multi-decadal periods for variables slowly changing in time, for instance the root-zone soil moisture (see bottom panel of Fig.3). Finally, it is important to emphasize 665 that an exhaustive evaluation of all land variables simulated in ERA5-Land is not feasible in a single paper. While this paper provides significant validation elements to demonstrate the added value of ERA5-Land, the wider scientific community is invited to carry out more detailed evaluation of individual components. For example, a more extensive validation of soil moisture following internationally agreed best practice  is highlighted as an area for further research.
While, in wider terms, we recommend the use of ERA5-Land fields over ERA5 for all type of land applications, one should 670 factor in their choice elements such as available computer and data handling resources, importance of spatial resolution versus data volume, area of application, temporal consistency, etc. With the public release of ERA5-Land, research and development has already shown some caveats of the dataset, such as the treatment of soil temperature in permafrost regions (Cao et al., 2020).
Further studies comparing with alternative well-referenced sites, as well as with regional and global datasets is encouraged.
For instance Pelosi et al. (2020) compared UERRA regional reanalysis (Copernicus Climate Change Service, 2020) forced by 675 ERA-Interim and ERA5-Land to assess the performance of evapotranspiration estimates based on weather data in the South of Italy.

Perspectives
The ERA5-Land dataset presented in this paper is the first operational land reanalysis of the ERA series. This paper has demonstrated its added value by comparing ERA5-Land estimates to a wide range of in situ observations, ERA-Interim and ERA5 680 reanalyses. Despite the overall observed improvement of the land states in ERA5-Land, there is scope for improvement for several components, which will be the focus for the construction of future enhanced versions. They are discussed below.
In the context of C3S, all operational products require an estimate of the associated uncertainty. Currently, estimates of uncertainty of ERA5-land variables are those corresponding to the ERA5 counterpart. ERA5 uncertainties are obtained by 685 running a 10-member EDA, which also provides the background-error estimates for the deterministic HRES 4D-VAR data assimilation system (Hersbach et al., 2020). First tests conducted at ECMWF running an ensemble of offline simulations with an ensemble of initial conditions and atmospheric forcing provided by the 10-member EDA of ERA5, indicate that the spread of the land surface variables is unrealistically low (not shown). The likely reason is that in the production of the ERA5 en-semble the physics of the surface model is not perturbed, i.e., the surface model is assumed to be perfect. To obtain realistic uncertainties, not only input and forcing parameters to the land surface model should be perturbed, but also key variables and parameters of the surface scheme (e.g. MacLeod et al., 2016;Orth et al., 2016), adding the contribution of land surface model error to the ensemble spread. Future studies are envisaged to investigate this path to provide meaningful uncertainties to land reanalysis.
One of the most important driving variables of ERA5-Land is precipitation. Systematic model-based precipitation biases can potentially spread into the land state estimates. ERA5 precipitation (used as input forcing of ERA5-land) benefits from a much improved data assimilation system assimilating millions of extra observations compared to its predecessor ERA-Interim.
Indirectly, ERA5-Land benefits from these extra observations as well. However, ERA5 still has large precipitation biases, especially in tropical regions. A recent bias-adjusted dataset based on ERA5 has been produced for impact studies with a 0.5 • resolution (WFDE5, Cucchi et al., 2020). That study demonstrated the added value of the bias corrections on large-scale hy-700 drological modeling. Future versions of ERA5-Land could consider similar approaches taking into account such bias-corrected coarser resolution forcing datasets as well as their availability in near-real-time. In addition to coarse-scale near-real-time biascorrection, high resolution downscaling could be also explored. Such a correction could be based on a climatological rescaling of precipitation based on a high resolution reference climatology (e.g. Karger et al., 2017) The land surface model used in ERA5-Land, CHTESSEL, also has the option to estimate carbon fluxes and its coupling with 705 plants transpiration through the A-gs formulation (Jacobs et al., 1996 ;Calvet et al., 1998;Boussetta et al., 2013a). This module is operationally active for the vegetation and it allows estimates for the carbon fluxes in a modular approach with evaporation being computed through a resistance approach (Jarvis et al., 1976). This choice is adopted since its integration in the operational numerical weather prediction still requires further developments (Boussetta et al., 2013a).  Chevallier et al. (2010) to correct for the biases in the 10-day budget of the modelled carbon fluxes at continental scales. This is currently being extended for usage in the FLUXCOM product (Jung et al., 2020) in order to bias correct the two components of the biogenic fluxes (Gross Primary Production and ecosystem respiration) separately. Figure  Another important research path with large potential is the revision and use of dynamic auxiliary data. Currently, ERA5-

720
Land assumes a fixed land cover, whereas Leaf Area Index (LAI) and albedo are based on a static monthly climatology. The former assumes that land cover remains unaltered for the complete reanalysis period and that cities are non existent, whereas the latter will not be able to accurately represent more frequent LAI anomalies. Based on the research conducted in the ESA-CCI programme, C3S provides (through the Climate Data Store (CDS)) climate data records of land cover at yearly frequency, as well as a close to real time global maps of LAI with 10 to 20 days latency. Moreover, recent studies have identified errors in the diurnal cycle of land surface temperature over Iberian Peninsula associated with the current land cover used in CHTESSEL and 1 hour model time steps (Johannsen et al., 2019;Nogueira et al., 2020). Therefore, the revision of the land cover and LAI to a new database as well as the introduction of their inter-annual variability is expected to provide more realistic land conditions as input to the surface scheme. There are also ongoing efforts to revise the vertical discretisation of the vegetation roots and soil layers distribution (e.g. Mueller-Quintino et al., 2016;Stevens et al., 2020), which could complement the land cover and 730 vegetation updates.
Finally, the coupling of the offline simulations to an offline data assimilation system is a promising approach. The advantages are multiple, in particular the assimilation of local datasets (precipitation radar observations, local soil moisture networks, etc.), permitting small adjustments of land state estimates. Undoubtedly, this will go accompanied by an increase in computational cost. Ongoing efforts to improve the parallelization of the land model will allow global high resolution simulations for several 735 decades to be performed at affordable computational costs.
All the components discussed above are currently under research, but they are not the only ones. ECMWF is currently working with a flexible, modular system called ECLand (Boussetta et al., 2021) which allows to develop separately several modelling aspects of the land surface, such as the increase of the number of soil layers or the introduction of a multi-layer 740 snow scheme, as well as progressing on other longer term perspectives such as the introduction of a groundwater storage or the reduction of the model time step. All the above ongoing developments will provide the basis for a future new version with improved accuracy for the land states at multi-decadal timescales.
7 Data availability ERA5-Land data are available through the C3S CDS, and at the time of writing this paper the data are available from January 745 1981. The data are accessible either via the user interface (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land? tab=overview) or through the CDS Application Program Interface (API). The data are updated with 2-3 months delay with respect to real time. However, a close to real time facility is planned to be implemented in 2021. The atmospheric forcing used to drive the land simulations is also available and already interpolated to the ERA5-Land grid. Note that in the CDS, and for user convenience, the data have been interpolated to a regular lat-lon grid of 0.1 • resolution.

750
All ERA5-land fields are made available at hourly temporal resolution and post-processed as monthly means. Monthly means are easier to handle and faster to retrieve, which is especially well suited for climate studies and also addresses an important requirement of reanalysis users. Two types of monthly means are post-processed: monthly means of daily means (average over all the hourly fields in a month) monthly means of synoptic means (averaged for a specific time of the day in a month) with < x > the estimate of the monthly average for the field x, N d the number of days in a month, N s the number of forecast steps in a day, M the forecast model, d the day and s the forecast step from 00UTC in a 24 h cycle. Note that the water and energy fluxes are accumulated from the beginning of the forecast time at 00UTC with a maximum of 24 h accumulation period.

760
For example, the monthly mean of surface runoff at 12UTC will provide the monthly averaged runoff accumulated from 00 to 12UTC. Therefore, for water and energy fluxes, Eq. 6 becomes: monthly means of daily means for water and energy fluxes Further technical details are also provided in the online documentation (https://confluence.ecmwf.int/display/CKB/ERA5-Land% 765 3A+data+documentation).
ERA-Interim surface data used in this study is freely available through the ECMWF catalogue: https://apps.ecmwf.int/datasets/ data/interim-full-daily/levtype=sfc/ ERA5 hourly data on single levels used in this study is also freely available through the C3S CDS (doi: 10.24381/cds.adbb2d47).

770
All soil moisture data used for validation is available through the International Soil Moisture Network, https://ismn.earth Access to the lake data used for evaluation in this paper is accessible as follows: -The Alqueva reservoir hourly data were provided by the Portuguese University of Evora, data are provided by demand, 897575 (see also Ménard et al., 2019). The snow depth dataset from the GHCN-daily network is available at https://www.ncdc.

790
GLEAM data were accessed from https://www.gleam.eu/ and the FLUXNET2015 Tier2 data set can be accessed from the FLUXNET data portal at https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/. This work used eddy-covariance data acquired and shared by the FLUXNET community, including these networks:  Author contributions. JMS produced the ERA5-Land dataset, designed the study and wrote the paper. ED, GB, SB, HH, DM, CB and JNT revised the paper and the figures, wrote some text and provided a critical review of the paper, CA pre-processed soil moisture data and NJRF conducted the corresponding analysis, GA processed and analysed the snow data, MC collected, processed and analysed all the lake data of this study, SH and EZ provided the evaluation of the river discharge data, BM processed and evaluated all the surface energy fluxes and run 810 the GLEAM model, MP collected and evaluated the skin temperature data, AAP processed and produced the bias-corrected carbon fluxes example. SH also proofread most of the paper.
Competing interests. To our knowledge there is no conflict of interest of the material presented in this paper.