ECLand: The ECMWF Land Surface Modelling System

Full text

Turn on search term navigation

1. Introduction

Land surface models (LSMs) describe the exchanges of the global energy, water and carbon fluxes with the land surface and improve our monitoring capacity to predict natural resources and their evolution in time. Besides providing boundary and initial conditions to Numerical Weather Prediction (NWP) systems and Earth System Models in general, land surface models also provide the terrestrial component for several research and operational applications depending on the represented processes and their degree of complexity. LSMs are nowadays used in local to global scales and sub-hourly to centennial scales through multiple applications going from climate change to hydrological and biogeophysical cycles such us floods and droughts or agriculture and CO $_{2}$ monitoring [1,2,3,4].

The ECMWF land surface scheme [5,6,7] is an essential component of the Integrated Forecasting System (IFS) but form also a self-standing model that can be integrated with specified meteorological forcing. In recent years, owing also to modularity of testing in offline and coupled applications, the land-surface modelling activity at ECMWF has been substantially and continuously expanded. An improved soil hydrology [8], a new snow scheme [9] and a multi-year satellite-based vegetation climatology [10] have been included in the operational IFS. These have had a positive impact on both the global water cycle and near-surface atmospheric variables compared to the original TESSEL (Tiled ECMWF Scheme for Surface Exchanges over Land) scheme used in ECMWF ERA-40 reanalysis [11] and in the ERA-Interim reanalysis [12], compared to HTESSEL used in ERA5 [13] and ERA5Land [14]. In particular the quality of seasonal predictions during extreme events associated with soil moisture-precipitation feedback was enhanced with the soil hydrology as in the European summer heat-wave in 2003 [15]. The new snow scheme resulted in a substantial reduction of near-surface temperature errors in snow-dominated areas by improving the thermal energy exchange at the surface (i.e., northern territories of Eurasia and Canada) [16]. More recently, the fixed maximum vegetation Leaf Area Index (LAI) field was replaced by a monthly climatology for LAI. This resulted in a reduction of near-surface temperature errors in the mid-latitude and tropical areas, particularly evident in Boreal spring and summer. The bare ground evaporation has been also enhanced over deserts and sparse vegetated areas by considering a lower stress limit in agreement with experimental findings (e.g., [17]) and results in a more realistic soil moisture for dry-lands. These model components are schematically illustrated in Figure 1 along with prospective developments.

Beyond the important day-to-day experience of running the land-surface scheme in an operational NWP environment, the participation in international projects such as AMMA (African Monsoon Multidisciplinary Analysis [18]) and GLACE2 (Global Land-Atmosphere Coupling Experiment-2), helped to better understand the mechanisms and areas of strong coupling between the land surface and the atmosphere based on experiments where the ECMWF IFS model was coupled with a realistic set of soil moisture fields.

With a continuous increase in spatial resolution and more broad applications of LSMs, as well as the growing availability of a range of diverse observations, a clear need for a more complex and accurate process representation is emerging. This continuous research effort aims at improving the land surface model parametrization for better surface prediction and ultimately better NWP and atmospheric composition forecasts. The current improvements are targeting different components of the land surface model, namely, the vegetation parametrization, the snow component, the soil hydrology, the open water/lake parametrization and also the representation of urban areas.

An efficient and scalable land surface system is essential for a routine NWP and service environment and can facilitate a paradigm shift for an Earth-Observation driven model development (Balsamo et al. [19]) under an Earth System framework, while linking with new possibilities brought by machine learning and Artificial Intelligence (AI) tools (Reichstein et al. [20]). The ECMWF land surface model has evolved as a modelling system (ECLand) that would facilitate such modular extensions for the benefit of efficient developments and external collaborations. Figure 1 illustrates the main ECLand system components and its prospective evolution elements that would facilitate its interaction with external projects and modular added services and being part of the full Integrated forecasting system. The ECland system attempts to resolve several processes linked to the water, energy and carbon cycles under an earth system modeling concept. Examples of theses processes are depicted in the different boxes of Figure 1 and detailed in the following sections. Important criteria for incorporating land-surface developments into ECLand are: Observability, for example, EO-data driven Calibration/Validation and Initialization, and added complexity if there is a route to successfully initialise and verify; Relevance for NWP and Copernicus applications, for example, leading to a substantial product enhancement; Impact on coupled Earth System forecast performance and key products, for example, improving NWP, extended range or seasonal atmospheric forecast skill, or providing skilful products for Copernicus services.

There is a long track record of participation in model intercomparison projects facilitated by WMO (e.g., GEWEX, WGNE), which above all verifies and validates different approaches, and sets standards for model development while identifying and publicising systematic model errors, model uncertainty and shortcomings. It also illustrates needs for added complexity and thus makes the case for any of the above criteria, namely observability, relevance, or impact. Meaningful examples are PLUMBER [21,22] dedicated to benchmarking, or the EartH2Observe FP7 project [23] where ECMWF provided high quality standardly formatted global meteorological forcing and coordinated a first multi-model Water Resources Reanalysis based on ERA-Interim/Land experience [24]. For hydrology, C3S assets such as ERA5L [14] motivated the ULYSSES project within the C3S sectoral hydrological applications with an offline multi-model stream. For Carbon, the CHE Horizon 2020 project [25] within the ECMWF CO $_{2}$ activities, have established ECMWF as a European CO $_{2}$ monitoring integrator and will within COCO $_{2}$ project coordinate a land multi-model benchmarking activity. This will be supported by experience gained within the Numerical Weather Prediction, for which a set of land-atmosphere coupling development priorities have been defined [26].

This paper describes the ECMWF land surface modelling system ECLand and a snapshot of its new developments. ECLand is based on the previous CHTESSEL parametrization with recent new developments focusing on the improvement of NWP within an earth system modeling concept. The recent developments aim at supporting new research questions in weather, hydrology, climate and environmental applications, ultimately at global kilometric-scale spatial resolution.

The configuration used in ECMWF’s operational NWP is described in the next section (Section 2: model description), followed by ongoing developments aimed at improving the model and its application under the model developments Section 3. An evaluation of the different new components is then performed, depending on the nature of the change with offline and/or coupled simulations, in Section 4 (evaluation of model developments) and a final Section 5 presents a summary and future perspectives.

2. Model Description

This section introduces the current land surface parameterisation operationally used at ECMWF. It describes the representation of the land surface fluxes of energy, water and carbon, and where appropriate, corresponding sub-surface quantities. The main surface scheme parameterisation structures were set with the Tiled ECMWF Scheme for Surface Exchanges over Land [5,6,7]. Each grid-box at the interface between the surface and the atmosphere is fractioned into tiles. There are currently seven tiles over land (high and low vegetation, bare ground, shaded and exposed snow, intercepted water and lake) and up to two tiles over sea and freshwater bodies (open and frozen water). Each tile is characterised with its own properties identifying separate water and heat fluxes and allow solving an energy balance equation for the tiles skin temperatures and taking into account physical processes that limit evaporation of vegetated areas and bareground (see the Supplementary Materials for detailed formulation).

One of the important characteristics of the ECMWF land surface scheme is the skin layer fast response, that allows for an implicit solving of the skin temperature equation in conjunction with the boundary layer’s vertical turbulent transport, as a natural integrated way within the IFS. Given the interdependence of the tiles’ skin temperature, the Best et al. [27] coupling approach is therefore used to jointly solve the skin temperatures equations in a fully implicit way.

The uniform global soil texture and the lack of surface runoff are shortcomings that have also motivated a revision of the land surface hydrology. In this context, Balsamo et al. [8] implemented a new runoff and infiltration parametrizations that are related to the standard deviation of orography and the soil texture. The snowpack is characterised by considering, snow density, thermal insulation properties interception of rain and the evolution of the albedo due to metamorphism aging processes [9,28].

The skin temperature over land is thermally in contact with a multi-layer soil (currently four layers, Section 3.3) or with a snow layer mantle overlying the soil, if snow is present (Section 3.2). The combined effect of melt energy, basal heat flux and top energy fluxes triggers the snow temperature variation. A modified Fourier diffusion equation that considers the soil water phase changes’ thermal effects is used to characterise the soil heat budget. Considering a zero-flux as the bottom boundary condition and a net ground heat flux at the top, allow then to solve the energy equation.

The snowmelt contributes to surface runoff and infiltration by depleting the snow mantle that collects the Snowfall. The rainfall is partitioned in fractions that are collected by an interception layer and the snow mantle, and a remaining throughfall fraction. This latter fraction contributes to surface runoff and infiltration. The multi-layer soil discretization which is shared for heat and water budget equations is used to solve Darcy’s law that govern the subsurface water fluxes. The bottom layer boundary assumes a free drainage condition and for the top boundary, surface evaporation and infiltration are assumed. While an additional root extraction sink from vegetated areas is considered in each layer. A revised bare soil evaporation enhances the soil-atmosphere water exchange especially over deserts and sparsely vegetated regions [29] and introducing a vegetation growth seasonal evolution allows to better modulate the evapotranspiration [10].

As an important component of the global carbon dioxide (CO $_{2}$ ) budget, a land surface CO $_{2}$ exchange module has been added to ECLand [30] enabling environmental forecasting applications which also involves interaction with atmospheric CO $_{2}$ concentration. The release of land biogenic CO $_{2}$ and the photosynthesis processes fixing carbon dioxide into biomass are parametrized in ECLand and respond to meteorological and climate forcing and to the natural biomes distribution including their stress conditions. The soil respiration is also parametrized in an NWP adapted way as a function of land-use. The photosynthesis and transpiration parametrizations are made modular within the land surface scheme allowing an independent interaction with the atmospheric CO $_{2}$ concentrations for global monitoring and prediction purposes.

Finally, a mixed layer model [31,32,33] is used to parametrize lakes and subgrid water as open water tiles. A specified surface temperature is assigned to open ocean, and in ocean coupled simulations, surface ocean temperature is driven by the ocean model tendency. A sea-ice cover specifies the fraction of the gridbox occupied by sea-ice, which is in thermal contact with underlying freezing ocean. The sea-ice temperature evolves using the heat budget equations of a four-layer ice model. While the lake model treats prognostically coastal and lakes’ ice temperature following a single layer ice scheme with an underlying water at freezing temperature.

The CO $_{2}$ -enabled land surface scheme, which includes the lake and subgrid water treatment, is the version currently used in the operational IFS system (CY47R1) and is also the base of the ERA5Land land reanalysis [14]. For the sake of completeness, a detailed mathematical description of its different components is provided in the Supplementary Materials and the related cited papers. In this paper, the model configuration based on operational IFS cycle (CY47R1) supporting daily forecasting from June 2020 will be referred to as the “current” model version.

ECLand encompasses the current model version, operational in 2021 and results from research developments that are foreseen entering in upcoming model upgrades and detailed in the next section.

3. Model Development

A continuous effort is made to improve the realism of surface-atmosphere interactions and to keep the model up-to-date with regard to state-of-the-art developments. The efforts described below are motivated by (i) the need to reduce biases and deal with compensating errors within the coupled NWP framework [26], (ii) the need to improve process-based parametrizations given the continuous increase in spatial resolution and the broader application spectrum, and (iii) the need to use more accurate and up-to-date ancillary data to deal with very high spatial resolution and reduce representativity errors. Hereafter, we introduce ongoing research work with different levels of maturity which are available in ECL and as configurations that can be activated independently or jointly. These developments are also categorised as ready for operational implementation or operationally viable for upcoming cycles.

3.1. Land Biosphere Representation

The vegetation developments that compose the land biosphere model are twofold, one focuses on the upgrade of the vegetation related maps, the other focuses on the improvement of the vegetation dynamics in time.

3.1.1. Land Use/Land Cover and Leaf Area Index

Current Land Use/Land Cover (LU/LC) are from the dated GLCCv1.2 data [34] which is based on observations from the Advanced Very High Resolution Radiometer (AVHRR) covering 1992–1993. An upgrade to more accurate and up-to-date vegetation LU/LC maps is introduced with a new automated framework that enables an easy use of the latest maps. These new maps are based on the ESA-CCI LU/LC brokered through the Copernicus Climate Change Services (C3S), which provides consistent maps at 300 m spatial resolution on an annual basis from 1992 to present (currently 2019). The maps have a total of 22 land cover classes based on the land cover classification system (LCCS) developed by the United Nations Food and Agriculture Organization (UN FAO). To be used in the model, the 22 LCCS classes are converted to the Biosphere-Atmosphere Transfer Scheme (BATS) classes through a new cross-walking table which allows removing “residual” (non pure) vegetation types such as the interrupted forest (which covers about 25% of the land points in the current ECMWF maps). It is important to note that having actual/pure plant types allows a better characterization of the model parameters, especially for parameters based on observed quantities.

The vegetation types in ECLand are categorised as low for crops, grass and herbaceous types and high for forest and trees types and their corresponding climatological cover is extracted form satellite based observations such as the the GLCC or ESA-CCI LU/LC products (please see the Supplementary Materials for the detailed equations and types classification). The difference in high and low vegetation covers between ESA-CCI and GLCCv1.2 based maps (Figure 2) shows an increase of low vegetation cover at the expense of the high vegetation cover in forest areas and a decrease of low vegetation cover favouring more bareground, especially over the Siberian and Australian deserts. The vegetation types show a split of the mixed and interrupted forest types between pure forest and low vegetation types. These differences suggest a substantial impact on the energy fluxes which will be explored in the evaluation section.

Additional developments are focused on improving the model phenology and its seasonality, which is represented in the model by the Leaf Area Index (LAI). To be used in the model with the tiles concept, the LAI needs to be split into low and high vegetation components. An update of this LAI disaggregation operator is developed. The current operator tends to overestimate the LAI during the transition periods (spring and autumn), while the new operator is more conservative of the observed LAI. The currently used operator relies on a LAI lookup table that represents the yearly maximum LAI per vegetation type and applies a seasonal cycle based on a MODIS LAI climatology [10]. The new operator extracts a seasonally varying LAI lookup table per vegetation type. In combination with the vegetation types map it defines the ratio of low and high LAI components which is then applied to disaggregate the satellite based climatological LAI. The revised operator allows thus to derive low and high LAI climatology components which when combined tend to conserve the observed total LAI climatology, in particular for the transition seasons. It also prepares for an operational assimilation of the LAI to allow better monitoring of the surface conditions inter-annual variability and representation of extremes as found in Boussetta et al. [35]. Following the changes in the high and low vegetation coverages, Figure 3 shows that for July the new LAI components tend to increase for high vegetation by up to 80%, while for low vegetation it mostly decreases by up to 60% in low vegetation dominated areas. In January (not shown) both low and high vegetation LAI tend to decrease by up to 80% compared to the current values.

3.1.2. Dynamic Vegetation

Taking into account vegetation dynamics in ECLand involves two different steps: in the first, dynamical vegetation cover is considered, and for the second, a prognostic LAI parameterisation is developed.

Seasonal Vegetation Cover

Unlike albedo and LAI, for which seasonality is taken into account, the current vegetation cover is implemented as a fixed map with no seasonal or inter annual variability following:

(1) $c = A C_{v e g} (v e g T y p e),$

where c is the total vegetation coverage, A is the climatological vegetation cover from satellite based land use/land cover maps and

C_{v e g}

is a vegetation type dependent density (see the Supplementary Materials for further details. In the new development, the seasonality of the vegetation cover is introduced by replacing the vegetation type dependent density and linking it to the LAI seasonality through the Lambert-Beer formulation (Figure 4), considering the vegetation clumping [36]. The vegetation density

C_{v e g}

is expressed then as:

(2) $C_{v e g} = 1 - e^{- k \cdot L A I},$

where k is a coefficient dependant on the vegetation clumping for low or high vegetation. Typical values of k range between

0.4

and

0.8

, here we consider a value of

0.6

, which was also adopted by Nogueira et al. [37]. In the latter study, when the clumping was also applied to the high vegetation, coupled simulations indicated a strong impact over the Northern Hemisphere circulation which require further investigation, therefore the implementation in this study is so far limited to low vegetation.

Prognostic Vegetation Parameterisation

The Leaf Area Index is simulated with a simple growth model based on the daily carbon balance at the leaves level following Gibelin et al. [38] and Calvet and Soussana [39]. The LAI is computed from the leaf biomass B derived from the net canopy assimilation $A_{nI}$ . Growth is the result of carbon accumulation obtained from assimilation of atmospheric CO $_{2}$ and senescence related to the photosynthesis deficit due to environmental factors. LAI is obtained from the product of the green leaf biomass B with the Specific Leaf Area ( $S L A$ ) which is a function of leaf nitrogen concentration dependent on vegetation types following:

(3) $L A I = B \cdot S L A = B \cdot e \cdot N_{l} + f,$

where

N_{l}

is the leaf Nitrogen concentration, and e, f are respectively the slope and intercept nitrogen plasticity parameters which vary by vegetation types. For global scale applications and to reduce drifts caused by forcing errors, the computed biomass is updated by a relaxation formulation toward climatological biomass values

B_{C}

(based on climatological LAI), the final biomass value

B_{f}

follows:

(4) $B_{f} = r_{B} \cdot B + (1 - r_{B}) \cdot B_{C},$

where

r_{B}

is the relaxation coefficient. For the sake of completeness, the growth model is further detailed in the Supplementary Materials.

3.2. Updated Snow Parametrization

A multi-layer snow scheme (ML) is introduced in ECLand, replacing the previous single-layer (SL) snow scheme. ML is an intermediate complexity snow scheme following the definition by [40], representing the vertical structure and prognostic evolution of snow temperature, density, liquid water content and albedo. A schematic of ML is shown in Figure 5. A detailed description of ML is reported in Arduini et al. [41] and in the following only the main characteristcs of ML are reported.

The heat budget of the snowpack is solved in two-steps. At first, a preliminary snow temperature is computed solving the heat budget in absence of melting, considering the vertical diffusion of heat and the absorption of solar radiation by each layer. For the latter, a parametrization adapted from Jordan [42] is used, in which the penetration of solar radiation decreases exponentially with depth and the extinction coefficient is a function of snow density. Then, melting and freezing processes are diagnosed by computing the heat content of each snow layer, also considering the liquid water content of each snow layer, and the snow temperature is modified accordingly to the latent heat released/absorbed. The movement of liquid water between layers is taken into account using a bucket-type formulation [16]. The heat flux at the snow-soil interface is computed as the sum of the residual heat flux at the end of the heat budget computation and the basal heat flux ( $G_{B}$ in Figure 5).

The number of active layers (N) varies depending on snow depth $D_{s n}$ , from a minimum number of 1 layer to a maximum of $N_{m a x} = 5$ . Multiple layers are used when the grid box is completely snow-covered, that is for $D_{s n} \geq 0.1$ m. As opposed to Arduini et al. [41], a different discretization algorithm is applied over flat terrain and complex terrain regions, the latter defined as regions with the standard deviation of the sub-grid-scale orography greater than 50 m. A variable vertical discretization is used to maintain a relatively high vertical resolution for the snow layers responding to fast time scales also for a deep snowpack. Over flat terrain, the thickness of the topmost layer in contact with the atmosphere is fixed to $0.05$ m; the second and third topmost layers, as well as the one in contact with the soil underneath, can increase up to a certain maximum $D_{i}^{m a x}$ as described in [41]; for $D_{s n} > 0.45$ m, the layer $N_{m a x} - 1$ is used as an accumulation layer. An example of the vertical discretization of a $1.0$ m-thick snowpack is shown in Figure 5b.

Over complex terrains, for $D_{s n} < 0.25$ m the same vertical discretization as described for a flat terrain is used. For a snowpack with $D_{s n} \geq 0.25$ m, the minimum ( $D_{m i n, i}$ ) and maximum ( $D_{m a x, i}$ ) thicknesses of each layer (i.e., the minimum and maximum allowed thickness for an active snow layer) varies with snow depth as follows:

$D_{m i n, i} = \min (0.25, 0.10 + α_{o} (D_{s n} - 0.25)) if i = 1 \dots, N_{m a x}$

$D_{m a x, i} = \{\begin{matrix} \min (0.25, 0.10 + α_{o} (D_{s n} - 0.25)) & if i = 1 \\ \min (0.30, 0.15 + α_{o} (D_{s n} - 0.25)) & if i = 2 \dots, N_{m a x}, \end{matrix}$

with

α_{o}

set to

0.10

. This implies that for a thick snowpack over complex terrain the snow layers are thicker than a snowpack with same depth over a flat region. A full description of the vertical discretization algorithm is detailed in Arduini et al. [41].

A number of snow physical processes are also changed in ML compared to SL. The snow heat conductivity is parametrized following Calonne et al. [43], with an additional term considering water vapor diffusion effects [44]. Changes in snow density due to wind blowing are an important process in particular at high latitude [45,46]. In ML a parametrization of wind-driven compaction is introduced, using the formulation proposed by Decharme et al. [46]. Snow cover fraction and albedo in ML are parametrized using the same formulations as in SL [9]. Modifications of these parametrizations will be the focus of future reaseach within ECLand.

3.3. Updated Soil Parametrization

With higher spatial resolution and more accurate grid parameters characterisation, in addition to broader applications and the need to improve the water budget in the model and coupled hydrological components, the importance of the soil depth and vertical discretisation within land surface models has increased. It affects the way land surfaces store and regulate water, energy and also carbon fluxes. The amount of water in the soil and its vertical distribution in the column are important for the regulation of heat and water vapour fluxes towards the atmosphere spanning a range of time scales from minutes to months in the coupled land-atmosphere system. This also depend on the soil parameters and is further modulated by the vertical roots distribution, and soil moisture stress function, which control evapotranspiration under soil moisture stress conditions. Currently in the ECMWF land Surface Scheme the soil column is represented by a fixed 4-layer configuration with a total of 2.89 m depth (7 cm, 21 cm, 72 cm, 189 cm). In the ECLand system we introduce new configurations with increased soil depth (up to 8 m with fixed layers depths) and higher vertical discretisation (up to 10-layers; 1 cm, 2 cm, 4 cm, 9 cm, 12 cm, 30 cm, 42 cm, 100 cm, 200 cm, 400 cm), including a dissociation between the treatment of water and heat fluxes (Figure 6) following [47]. In these configurations, the upper layers which have a thinner depth, allows a better match with satellite microwave soil moisture observations.

An optional formulation for the distribution of roots density has been also developed in ECLand. It is based on a uniform distribution, which is appealing due to its simplicity and straightforward first-order assumption. This uniform distribution replaces the current exponential formulation and depends only on the maximum rooting depth ( $R_{D M A X}$ ). This root distribution formulation was investigated by Stevens et al. [48] and showed potential for improvement of the water and energy fluxes. The transformation of root zone soil moisture to water stress is also an open question, and there are different approaches to relate root zone soil moisture and the water stress factor [49,50]. In ECLand it is also possible to change the linear transformation between root zone soil moisture and water stress by a non-linear relation, which showed a large sensitivity in offline and coupled forecasts [51]. Such approaches, combined with the rooting depth distribution, are likely to further enhance the model capability to represent root water uptake, but require further evaluation and calibration. The quality of the soil ancillary data and its related parameters is also essential for the regulation of water, energy and carbon fluxes. The coupling of ECLand with the multiscale parameter regionalization system (MPR [52]) would facilitate the integration and calibration of such parameters related to new soil maps.

3.4. Surface Water Representation

An upgraded versions of the physiographic datasets used in ECLand to generate the land sea mask and the lake parameters is described. In addition to the ecosystem data, the new land-water mask uses a newly generated glacier data and a new land-water and lake fraction mask. The lakes are also benefiting from updated lake depth data. The main data sources for land-water mask are from the Joint Research Centre (JRC) with nominal resolution 30 m (1 $^{″}$ ). The high resolution information at 30 m is used in the data merge then aggregated to 1 km.

The upgraded inland water mask is based on the upgraded fractional land-water mask at 1 km resolution. First, fractional land-water mask is translated into “land” and “water” binary mask, then by using an innovative “waterbody” separation algorithm [53] “water” is separated into ocean and inland water bodies. The algorithm is pixel-by-pixel and iterative based on the spatial window width, and the number of iterations which depends on the resolution and complexity of the coastline. At the end of all iterations a visual check is needed to correct freshwater coastal lagoons where needed. The Azov sea and Maracaibo lake are also processed as inland water given that they are better represented in the model by the lake parametrization which illustrates some of the assumptions required in this process. Figure 7 illustrates the main differences from the new land water mask and the current operational one (CY47R1). It shows that several inland water (lake) appear in the new map especially in north Canada and Russia while some lakes are removed from central Australia. Detailed description of the data sources and the processing chain for the updated surface water is provided in the Supplementary Materials.

3.5. Coupling with a Hydrodynamic Model

River discharge is a key component of the global water cycle, integrating meteorological and hydrological conditions on multiple temporal and spatial scales. Weather forecast uptake in hydrological forecasting [54] is a growing area, and developments by the climate modeling community to enhance the representation of the Earth system have motivated a tighter integration of atmospheric models and the surface water balance [55]. River discharge provides a spatially integrated measure of the surface water budget and can be used to evaluate forecast skill [56]. Furthermore, in atmosphere-ocean coupled numerical weather prediction, freshwater inflow to the ocean can impact ocean circulation Zuo et al. [57]. Therefore, the representation of rivers and their role on the ocean and land water is relevant in several applications, motivating the integration of the Catchment-based Macro-scale Floodplain model (CaMa-Flood) into the ECLand system.

CaMa-Flood model was developed to simulate continental-scale river hydrodynamics. The global river network (excluding Antarctica) is discretized in irregular unit-catchments, representing hydrological units. Sub-grid topographic parameters describing the river channel and floodplains profile are used to diagnose flooded area and water levels from water storage. The irregular-shaped unit-catchments are assigned to a rectangular grid, providing a “grid-vector hybrid river network map”. This approach allows a realistic representation of sub-grid topography while keeping data handling simple. River flow velocity is computed by the local inertial equation [58]. There is also an option to compute flow velocity using a kinematic wave approximation. The only prognostic variable, water storage, evolves in time according to the input runoff and upstream inflow and downstream outflow at each unit-catchment. The detailed description of CaMa-Flood can be found in the description papers [59,60,61].

The key characteristic of CaMa-Flood is the explicit representation of flood stage, representing water level and inundation extent, which complements the conventional use of river discharge [62]. This also allows exploring two-way coupling with the land surface model [63]. The baseline topography is based on MERIT DEM/MERIT Hydro [64,65]. River channel width is objectively parametrized using satellite measurement [66,67].

The CaMa-Flood source code is integrated in the ECLand system allowing stand-alone CaMa-Flood simulations or simulations with 1-way and 2-way coupling between ECLand and CaMa-Flood. The coupling is performed sequentially and the code is integrated as a single executable (without external couplers) following a similar approach as the coupling with the ocean model to IFS [68]. In 1-way coupling, the default configuration, surface and sub-surface runoff are sent to CaMa-Flood. In the stand-alone mode CaMa-Flood reads from input files the runoff data, reproducing the 1-way coupling results. This configuration has a lower computational cost, allowing for testing of the hydrodynamic model as well as driving the model with other runoff data sources like reanalysis. The 2-way coupling between ECLand and CaMa-Flood was also experimentally developed allowing initial testing of the interaction between inland water in ECLand and flooded areas in CaMa-Flood. Inland water evaporation (in ECLand) can be extracted from flooded areas in CaMa-Flood, and the flooded area extent can update the inland water fraction in ECLand. CaMa-Flood integration in the ECLand system is only available in offline mode supporting shared memory parallelism via OpenMP. A full integration of CaMa-Flood in the IFS is envisaged in upcoming cycles to support coupled atmosphere-ocean-land-rivers simulations aiming at a full representation of the water cycle. This will still require developments both at the level of distributed memory parallelism in CaMa-Flood and in the coupling of river discharge with the ocean model.

The river network description used in CaMa-Flood is currently available in several global regular latitude/longitude resolutions: $1 / 4^{\circ}$ , $1 / 10^{\circ}$ , $1 / 12^{\circ}$ , $1 / 20^{\circ}$ , and $1 / 60^{\circ}$ . However, each unit-catchment in CaMa-Flood has an irregular shape and ECLand can be setup on any regular or the IFS model’s Gaussian reduced grid. To accommodate the different spatial grids between these two components, the fluxes are interpolated using pre-computed weights that account for the area of interception of each grid element. For the 2-way coupling, the inverse of the weights is used to guarantee mass conservation. These weights are computed as part of the pre-processing of the system, which also allows the definition of regional domains, which are relevant for testing and for applications that do not require the global setup.

3.6. Urban Parametrization

A single-layer-type urban model is implemented within ECLand. A detailed description of the scheme, which includes a basic assumption of urban geometry, computes roughness lengths and considers fluxes from multiple surfaces (Figure 8), can be found in [69]. The Urban scheme is applied as a new sub-grid surface tile and uses mapped parameters and derived variables to compute surface energy and water fluxes. The urban environment is represented as both a canyon, consisting of walls and road, and a roof fraction. An infinite length canyon assumption is made to factor shadowing from buildings and is dependent on both urban parameters and the solar zenith angle. Building height and canyon characteristics are used to estimate surface roughness and any concept of canyon orientation is removed. Energy and moisture exchange between the surface and both the atmosphere and sub-surface is defined by urban parameters, including a simple representation of urban drainage. The scheme consists of both time varying (e.g., albedo based on solar zenith angle) and fixed variables (e.g., thermal heat capacity of a building). The variables required for the urban scheme are based on previous similar schemes [70,71,72] and listed in the Supplementary Materials.

Currently, ECLand consists of eight sub-grid surface tiles and a lake tile. To represent the urban environment we have introduced a tenth tile and updated the high vegetation snow to include an urban representation. When considering tile fractions over land, the value summed across all tiles is equal to 1. For the land surface this was previously a combination of wet land, low vegetation, high vegetation and bare soil and, for the latter three with and without snow cover. Over urban areas the addition of an urban land cover map would increase the sum above one. To prevent this, existing tile fractions are re-weighted to account for the urban tile following an update of the ECLand tiles Equation (5) where high vegetation ( $C_{H}$ ), low vegetation ( $C_{L}$ ), bare soil ( $C_{B}$ ) and snow under high-vegetation ( $C_{H_{snow}}$ ) covered fractions are updated from the original values ( $C_{H_{orig}}$ , $C_{L_{orig}}$ ) accounting for the urban fraction ( $C_{U}$ ).

(5) $\begin{matrix} C_{H} = C_{H_{orig}} \times (1 - C_{U}), \\ C_{L} = C_{L_{orig}} \times (1 - C_{U}), \\ C_{B} = (1 - C_{H} - C_{L} - C_{U}), \\ C_{H_{snow}} = C_{sn} \times (C_{H} + C_{U}), \\ C_{U_{dry}} = (1 - C_{wet}) \times (1 - C_{sn}) \times C_{U}, \end{matrix}$

where

C_{sn}

C_{U_{dry}}

and

C_{wet}

represent the snow covered, dry urban and wet fraction, respectively.

For land surface skin, future developments will merge the urban cover product with the remaining land cover so this re-weighting process is no longer required. For certain elements, for example, hydrology, this prevents a weighted average that neither fits an urban or non-urban regime. Currently, for sub-surface thermodynamic and hydrological properties, the urban tile re-weights the uppermost soil layer. The urban land cover used here is based on grid cell specific information provided by ECOCLIMAP-SG, at an initial 300 m horizontal resolution [73].

3.7. Model Efficiency

The surface model inherits the same “array-block” data structures as all other IFS physics parameterisation packages and is therefore capable of benefiting from modern processors and GPUs. These data structures allow for the exploitation of shared memory parallelism via OpenMP as multiple threads can work independently on distinct array-blocks, therefore leading to a reduction in memory consumption and improved load balancing, with the scaling of these blocks tailored depending on the target architecture. This makes ECLand a very well adapted and efficiently integrated component of the IFS. The offline ECLand model version equally supports distributed memory parallelism via the Message Passing Interface (MPI). The new implementation allows to significantly accelerate the computational performance of ECLand in offline mode at all resolutions and enables efficient simulation at very high resolutions such as 1 km (with a throughput of one year of global simulation per day). This throughput is important for a range of new developments including accounting for urban areas (see later Figure 26 in Section 3.6) and their carbon emissions within the future Copernicus CO $_{2}$ service as well as the foreseen developments under the Destination Earth initiative (DestinE), which aims at creating a high-precision digital twins of the Earth [74].

One of the largest scalability bottlenecks when running ECLand in offline mode at very high resolutions (e.g., 1 km) is input/output (I/O). Ongoing development would attempt to mitigate this by exploiting the capabilities of modern parallel file systems and higher level libraries such as NetCDF4. Parallel efforts focus to consistently adapt ECLand to ECMWF’s multi-IO development including optional output in different formats such as NETCDF or GRIB2.

Additional development towards increased model efficiency is the ability to run on accelerated hardware such as Graphics Processing Units (GPUs). An effort towards porting the surface model to run efficiently on GPUs has already started at ECMWF under the umbrella of the Scalability programme [75].

4. Evaluation of Model Developments

4.1. Data and Methods

New developments are usually evaluated with either offline or both offline and coupled forecasts depending on their target applicability and level of testing. For specific applications, only offline land surface modeling is needed while for others, coupled simulations are necessary, e.g., if the land-atmosphere interactions affect the simulation results. For implementation within the ECMWF operational IFS, a more strict assessment is adopted: new developments are evaluated offline, in coupled forecasts and data assimilation configurations. More recently, under the concept of earth system modelling, hydrology related scores, which introduce additional constraints for uncertain parameters, started to be also considered aiming at improving the core objective of NWP as well as other earth system applications. To characterize the impact of the new developments mentioned in this paper, the model performance is evaluated independently according to each developments’ level of maturity following the configurations described in Table 1. In this paper, hydrological verification is only performed for the coupling with the river discharge development. All 2D and 3D simulations are run at the Tco399 Gaussian octahedral grid (25 km), the offline simulations are driven by the ERA5 atmospheric forcing while the coupled forecast simulations are initialized with corresponding offline simulations to avoid unstable surface conditions related to spinup problems.

4.2. Observational Data to Evaluate Model Developments

The new developments described above are evaluated against the following observations or observation-based products.

4.2.1. FluxCom

The FluxCom dataset (www.fluxcom.org/, 1 February 2021) is used to evaluate the impact of the new developments on the energy and carbon fluxes. This dataset uses machine learning techniques, namely Random forests (RF), Artificial Neural Networks (ANN) and Multivariate Adaptive Regression Splines (MARS), to upscale FLUXNET in-situ flux observations in space and time and integrate them with remote-sensing and meteorological forcing data. It uses the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data for the FluxCom-RS version and additionally uses reanalysis meteorological forcing in the FluxCom-RS + METEO version. FluxCom provides fluxes for the RS + METEO version at 0.5° spatial resolution and daily temporal resolution for the 1979–2013 period and the RS version is available at 8-day temporal resolution with 0.0833° spatial resolution for the 2001–2015 period. Further details of the FluxCom data is available in Jung et al. [76].

4.2.2. GLEAM

The Global Land Evaporation Amsterdam Model version 3 (GLEAM v3b; Miralles et al. [77], Martens et al. [78]) is used to evaluate the impact of the vegetation updates on evapotranspiration. GLEAM is a semi-empirical-process-based model that computes evapotranspiration and its components globally over land at daily scale using a set of algorithms driven by satellite-based observations. Potential evaporation is derived over four land cover tiles (low vegetation, tall vegetation, bare soil and open water) using the Priestley-Taylor formulation. Then actual transpiration is calculated by multiplying the potential evaporation by a stress function based on root-zone soil moisture and vegetation optical depth. Soil moisture is the result of a soil-water balance model that computes the evaporative flux and is driven by observation-based precipitation. The last component of GLEAM is the interception loss which follows the Gash analytical model driven by precipitation and vegetation properties. Further details of GLEAM are documented in Martens et al. [78] and Miralles et al. [77].

4.2.3. FLUXNET

LaThuile FLUXNET data ( [79]) used are observations from 253 research sites (https://fluxnet.org/data/la-thuile-dataset/, 1 February 2021) which belongs to the FLUXNET eddy covariance network. Net ecosystem exchange (NEE), sensible heat flux (H) and latent heat flux (LE) at high temporal resolution (30 to 60 min) are provided in this dataset. Additionally derived gross primary production (GPP) and ecosystem respiration (Reco) are also provided based on the partitioning method proposed by Reichstein et al. [80]. Only data without gap filling and flagged as high-quality are used.

4.2.4. WOFOST

WOFOST is a crop growth model to estimate yields quantitatively. It is developed by the Department of Theoretical Production Ecology (Wageningen Agricultural University, The Netherlands) and the Center for Agrobiological Research and Soil Fertility (Wageningen, The Netherlands). This model is implemented currently within the MARS Crop Growth Monitoring System (CGMS), allowing the estimation of biophysical variables related with crop yields such as potential biomass production, crop development stage, and so forth. WOFOST raw simulations are done at soil unit level, that is, for each of the soil types within a soil polygon. These data were provided by JRC under the EU FP7 ImagineS project framework.

4.2.5. LSA SAF LST

Land Surface Temperature (LST) is taken from the Satellite Application Facility on Land Surface Analysis (LSA-SAF). LST is derived from measurements by the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board the Meteosat Second Generation (MSG). The LST retrievals are based on the generalized split-window technique [81] and are available every 15 min for all clear-sky pixels within the MSG disk with a resolution of 3 km at the nadir. Following previous work on model evaluation with satellite retrievals [37,82], the data is aggregated to a coarser grid (0.25 $^{\circ}$ × 0.25 $^{\circ}$ ) only keeping the valid pixels (fraction of clear sky). For the comparison with model simulations, clear sky conditions, synchronous in time are considered: (i) at least 70% of the fraction of valid pixels and (ii) total cloud cover lower than 30% [82]. For each day and pixel, when at least 75% of the 3-hourly data was available (6 out of 8), the daily maximum and minimum LST are extracted to compute the diurnal range.

4.2.6. ESA-CCI Soil Moisture

The ESA-CCI soil moisture product (https://www.esa-soilmoisture-cci.org/, 1 February 2021) is also used to evaluate the impact of the new developments on soil moisture. The ESA CCI soil moisture product is provided by the European Space Agency for the 1981–2019 period at 0.25 $^{\circ}$ spatial resolution and daily temporal resolution. It is generated with soil moisture observations from active (ERS1-2 SCAT and MetOp ASCAT A-B) and passive (SMMR, SSM/I, TMI, AMSR-E, WindSat, AMSR2 and SMOS) microwave satellite sensors. The combined product is obtained through the assimilation process of both C-band scatterometers active observations that use the change detection algorithm and the passive microwave observed land surface brightness temperature. Further details on ESA-CCI SM data and its generation procedure are found in Gruber et al. [83], Dorigo et al. [84] and Gruber et al. [85].

4.2.7. ESM-SnowMIP

Observations of snow depth, mass and density collected in the Earth System Modelling Snow Intercomparison Project (ESM-SnowMIP [86]) are important for the evaluation of performance of the new snow scheme at selected sites. ESM-SnowMIP consists of ten sites at which detailed snow-related observations are available for more than 7 years (some sites more than 15 years) continuously.

4.2.8. SYNOP Observations

SYNOP observations of 2-m temperature and snow depth play an important role in the assessment of the impact of the new snow scheme as well as the urban model. SYNOP observations are exchanged through the WMO Global Telecommunication Service (GTS) and are used everyday for initialisation and monitoring of NWP models. The usage of in-situ SYNOP observations within the ECMWF IFS model is described in Haiden et al. [87].

4.3. Evaluation of Updated Land Cover and Vegetation

4.3.1. Land Use/Land Cover and Leaf Area Index Evaluation

Evaluation of the new ESA-CCI LU/LC map and LAI disaggregation operator is performed through global offline and forecast coupled simulations. The offline simulations are driven by ERA5 meteorological data. A control simulation (CTR) with the current operational maps is compared with a simulation that includes the new ESA-CCI LU/LC map and the new LAI disaggregation operator (CLAI). The impact of the changes in LU/LC and LAI are assessed by evaluating the mean diurnal range of Land Surface Temperature (LST). The simulations are compared with the (LSA-SAF) LST. For each day and pixel, when at least 75% of the 3-hourly data was available the daily maximum and minimum LST are extracted to compute the diurnal range (representing the amplitude of the diurnal cycle) and the evaluation is performed for different seasons for the period 2004–2015. The mean LST diurnal range bias of the control simulation during June-August (JJA) and September-November (SON) for the 2004–2015 period is shown in Figure 9a,b. These two seasons were selected due to the high coverage of clear-sky conditions in the MSG disk, but the main features of the biases are similar in the remaining seasons. In both seasons there is wide-spread and substantial under-estimation of the diurnal range, which can reach values larger than 6 K in several regions such as the Iberian Peninsula, Southern Africa and East South America. During JJA there are indications of an over-estimation of the diurnal range in West Sahara and the Arabian Península (Figure 9a). However, in these two regions there are some concerns on the quality of the LST satellite retrieval due to contamination by aerosols (West Sahara) and emissivity directional effects (Arabian Península). The change in mean absolute error of the mean LST diurnal range between the simulation with the new ESA-CCI LU/LC and LAI (CLAI) and the control simulation for JJA and SON are shown in panels c,d of Figure 9. The results show a general reduction of the under-estimation of the diurnal range in the CLAI simulation. The increased diurnal range in CLAI is mostly driven by higher daytime maximum LST (not shown) and are consistent with previous results that indicated the added value of the updated LC and LAI in CLAI in terms of the LST diurnal cycle. Previous studies have also explored the role of vegetation seasonality [37,88], suggesting that further improvements of the LST diurnal cycle can be achieved.

In addition to the evaluation of the LST diurnal range, the mean evaporation of the CTR and CLAI simulations are compared with two independent estimates of evapotranspiration (ET) mostly constrained by remote sensing data: (i) the GLEAM version 3b (Martens et al. [78]) and (ii) FluxCom Remote sensing ensemble flux products [76]. Figure 10 shows the comparison of evaporation (ET) 2003–2015 annual mean between the control run (CTR) and the CLAI run. Left panels show the CTR percentage bias ( $p b i a s = 100 (e x p - o b s) / o b s$ ) with regards to GLEAM (GL3b) and FluxCom (FCRS) ET products as observation proxies (a and c panels respectively). The right panels show the absolute $p b i a s$ difference between the CLAI experiment and the CTR experiment, where GLEAM has been used for verification in panel (b) and FluxCom in panel (d). These differences show an overall reduction of the bias compared to both products (green colour). However, the results also show that on an annual time-scale, the improvements in $p b i a s$ with CLAI with respect to CTR are fairly small compared to the differences of the observational-based products. The latter is illustrated by comparing GLEAM and FluxCom in panels (a) and (c). It highlights the uncertainty in estimating ET as also confirmed in Jung et al. [76].

Evaluation of coupled forecast simulations is performed against the operational analysis. Given the longer time scales of the land surface processes compared to the atmosphere, and to overcome spin-up issues caused by unstable surface state, the surface initial conditions of each forecast run are obtained from the corresponding surface offline simulations. The synoptic situation is in the predictable range when focusing on short-range forecasts and an impact on near-surface atmosphere through sensible and latent heat flux and radiative forcing is also triggered by changes in land cover and LAI.

Screen-level temperature and dewpoint can be verified using the analysis, which draws closely to the surface synoptic observations. Figure 11 shows the 2 m temperature normalised RMSE difference between forecasts based on the new ESA-CCI LU/LC map with LAI disaggregation operator and forecasts based on the current operational maps with regards to the operational analysis for the 4 seasons.

A mixed signal is seen for different areas with a marked reduction of the 2 m temperature rmse (bluish colour) when using the updated vegetation especially in the northern hemisphere during MAM and JJA (panels (b,c) respectively), while JF and SON (panels (a,d) respectively) are marked with a mainly neutral impact. Although the overall results of introducing these updated vegetation maps are neutral to positive, a few areas still show a degradation in some atmospheric parameters pointing to the need for additional investigations and a parameter optimization with respect to the complex boundary layer turbulence interactions. This development is crucial for an accurate land surface representation and is therefore part of the operationally viable research for upcoming cycles.

4.3.2. Evaluation of Dynamic Vegetation

Seasonal Vegetation Cover Evaluation

The vegetation cover seasonality is also evaluated with global offline and forecast coupled simulations. Figure 4b shows low vegetation cover difference between January and July when using a clumping based seasonality. This is in contrast with the current static vegetation cover field, which has the same vegetation cover all year round. Figure 4a shows the difference of the seasonal bareground cover for July as compared with the current static map. This variability in cover impacts the global circulation as it is directly linked with the surface roughness. Moreover it affects surface albedo especially in winter and spring during the snow season for example, through the shading of snow under high vegetation with a strong impact on energy fluxes. This configuration is also assessed in Nogueira et al. [88] and showed an overall bias reduction of the maximum LSA-SAF LST with some scattered deterioration, especially in south America which was attributed to a misrepresentation of land use related parameters. In the forecast experiments, this change is tested on top of the previous ESA-CCI + LAI (CLAI) configuration. The impact on the near surface parameters in this case also showed mixed patterns which tend to amplify both positive and negative signals found when using the CLAI configuration only. Figure 12 illustrates the 2 m temperature normalised RMSE difference between the coupled forecast simulations based on the seasonal vegetation cover (FSCV) and the simulation based on the current operational static cover map (CTR). The operational analysis for the 4 seasons is used as a proxy for observations. The mixed signals in the verification emphasize that vegetation changes are not independent of other model settings which often require further optimisation of vegetation related parameters such us roughness length, vegetation minimum stomatal resistance and the surface skin conductivity. They all modulate the exchanges between the surface and the atmosphere and point to the need to further adjust this development as operationally viable for future cycles.

Prognostic Vegetation Phenology Evaluation

The prognostic phenology evaluation is performed with offline in-situ and offline global simulations driven by ERA-Interim forcing. A relaxation coefficient $r_{B}$ (Equation (4)) of 0.05 day $^{- 1}$ was used in this case. The results show that this simple growth model is able to detect growth anomalies with regards to a multi-year climatology. Figure 13 illustrates the above ground anomaly index ( $A I_{A G B}$ ) for November 2010 with regards to the 1999–2013 climatology defined as

(6) $A I_{A G B} = \frac{A G B_{N o v 2010} - A G B_{C L I M_{1999 - 2013}}}{A G B_{C L I M_{1999 - 2013}}} .$

November 2010 is chosen because it is characterized by a drought in the Horn of Africa and a wet event corresponding to drought recovery in central and eastern Australia. Given the tight link assumed in this parametrization between the above ground biomass and LAI, this map correlates well with both Copernicus global land (CGLS) observation-based and simulated LAI anomaly maps of the same period (not shown). In the Horn of Africa, the above-ground biomass can be less than 20% of the 15 year mean ( $A I_{A G B}$ < −80%), while for the Australia wet case, the anomaly is positive and in some areas it exceeds the 15 year mean by 60%.

The in-situ evaluation is performed with the Flux tower network (FLUXNET) data (Baldocchi et al. [79]) and the WOrld FOod Studies (WOFOST) products which are currently developed by the Crop Growth Monitoring System (CGMS) within the Monitoring Agricultural Resources unit (MARS) of the European Joint Research Centre (JRC).

Figure 14 is an illustrative example of this evaluation, where Figure 14a shows that the growth model reasonably simulates the LAI (black line) as compared to both the observed LAI (near real time Copernicus Global land LAI product GEOV1) and its corresponding 1999–2013 climatology (purple line). Although it slightly underestimates its maximum value for the considered Hungarian site (hu-mat), the correlation with the near real time LAI is slightly improved.

In this case, for 2004 at the hu-mat site, the gross primary production (GPP) based on prognostic LAI correlates better ( $r = 0.925$ ) with smaller RMSE ( $R M S E = 1.105$ ) than when using the GEOV1 LAI ( $r = 0.894, R M S E = 1.260$ ). The time series of the above-ground biomass (AgB) from 1999 to 2013 for the Hu-mat site is illustrated in Figure 14b. It shows that using the growth model (black line), the AgB correlates much better with the WOFOST products than when using the prescribed LAI climatology. These results are indicative of the overall impact of the growth model at the site level, although for other sites (not shown) the impacts may fluctuate. It should be noted here that this development is not implemented in coupled forecast mode and therefore still requires further testing to be included in the operational system.

4.4. Evaluation of Updated Snow

The new multi-layer snow scheme has been evaluated in offline simulations forced with observed meteorological data from the ESM-SnowMIP network [89], as well as in global offline simulations driven with ERA5 forcing data and in coupled ocean-land-atmosphere forecasts.

4.4.1. Site Evaluation at ESM-SnowMIP

Figure 15 shows summary statistics for Snow Water equivalent (SWE) and snow depth (SD) using the SL and ML5 snow schemes, evaluated at the ESM-SnowMIP sites. ML5 improves the bias, correlation and standard deviation of the errors at most of the sites considered. The largest improvements are found at Col de Porte (cdp), as reported already in Arduini et al. [41]. Senator Beck (snb) is the only site showing a consistent degradation in both SWE and SD, for which both correlations and errors increase for ML5, see [41] and discussion therein. SL has an overall positive bias in SWE and SD at most of the ESM-SnowMIP sites, at some sites exceeding 20%. ML5 generally reduces the biases, changing sign (negative) at some sites.

The impact of the new vertical discretisation on mountain sites is shown in Figure 16, by comparing a simulation at the Swamp Angel site with and without the snow depth-dependent vertical discretisation over complex terrains. The latter configuration, called ML5-noOro, corresponds to a simulation with the same vertical discretization of the snowpack as used by Arduini et al. [41]. ML5-noOro shows an overestimation of melting during the end of the snow season in particular years, which is alleviated in ML5. Considering the entire time period, the impact over mountain sites of the new parametrization is generally positive (not shown).

4.4.2. Global Offline Evaluation of the Snow Depth

The synoptic observations of snow depth on the Global Telecommunication System (GTS), which are routinely assimilated in the ECMWF snow analysis [90], are used to evaluate global offline simulations with the multi-layer (ML5) and single-layer (SL) snow schemes. The evaluation period covers the months between fall to spring (September to May) from 2014 to 2018. Figure 17 shows a general reduction of the root-mean-square error (RMSE) of snow depth when using the multi-layer snow scheme, with the largest improvements over North America and complex terrain regions. Analysis of the error distribution shows a reduction of the median value of the bias and RMSE, and generally an improvement at stations characterized by the largest errors (see Figure 17b).

4.4.3. Global Coupled Evaluation

Evaluation of the multi-layer snow scheme in coupled forecasts has been performed over multiple winter periods looking at the impact on weather parameters and surface fluxes using a wide range of diagnostics [41,91]. Using a thinner snow layer with a low thermal inertia enables better representation of strong cooling events occurring in snow-covered conditions, because the snow surface responds on quicker time-scales to changes in the atmospheric forcing. This leads to an increase of the amplitude of the diurnal cycle of 2 m temperature over snow-covered surfaces, as well as to an increase of its variability.

Figure 18 shows the absolute mean bias and root-mean-square error differences over a winter period (DJF) between coupled forecasts at a lead time of 2 days using the multi-layer snow scheme and the current single-layer snow scheme. Two metre temperature observations from the synoptic station network, used routinely in the ECMWF assimilation system for monitoring and diagnostics, are used as a benchmark.

Results indicate a reduction of the mean error in particular in the high-latitude regions, which can be associated to an improved prediction of very cold temperatures over snow-covered regions. The RMSE of T2 m is reduced in the snow multi-layer forecasts over most of Eurasia regions, however it increases over certain regions, for instance in Alaska and West Scandinavia. This can be partly associated to errors in other variables and processes, like cloud cover and/or cloud microphysics, for which the new scheme is more exposed more “responsive”, as described in [41]. However, the increase of RMSE over complex terrain can be also related to an increase of the preexisting biases in the control forecasts, which can be partly due to the relatively low horizontal resolution used in these experiments (25 km). The increase in variability can be advantageous in a probabilistic framework, through the increase in forecast spread [92]. The overall results of the updated snow scheme indicate its readiness for implementation in the next operational cycle.

4.5. Evaluation of Updated Soil Parametrization

The updated soil depth and vertical discretization is assessed with surface offline simulations forced with ERA5. A control simulation with the default 4 layers and 2.89 m depth is compared with 2 additional configurations; one using a 10-layer/8 m depth configuration for both water and temperature (OSML10_TW) and in the other, the depth is only extended for temperature while keeping the default 2.89 m depth for water (OSML10_T). Evaluated against the FluxCOM energy flux products, both 10-layer configurations show different signals for the latent and sensible heat fluxes compared to the control standard 4-layers configuration. As an illustrative example, Figure 19 shows that using the 10-layer configuration results in an increase of the sensible heat flux correlation with the FluxCOM product compared to the standard 4-layer configuration, while the latent heat flux shows a general decrease of the correlation with increase over arid areas. Note that in this case, the 10-layer configurations (OSML10_TW) and (OSML10_T) have very similar correlations for flux scores, although the configurations with 8 m for temperature and 2.89 m depth for water has different patterns in subsurface runoff (figures not shown).

Following the evaluation of the Land Cover and LAI impact on the diurnal range of LST (see Figure 9, the impact of the increased vertical resolution in the soil is also assessed by evaluating the LST diurnal range. Figure 20 shows the absolute bias difference of the mean diurnal range of LST between the OSML10_TW simulation and the control simulation for JJA (top panel) and SON (bottom panel). Similarly to the reduction of the diurnal range underestimation seen in CLAI (Figure 9), the increased vertical resolution of the soil near the surface also reduces the diurnal range under-estimation over the Iberian Peninsula, South Africa and Eastern South America. In general, the increased vertical resolution drives an increase of the LST diurnal range, due to a reduction of the thermal inertia of the top soil layer. While this compares favourably in most of the MSG disk, in the Western Sahara and Arabian peninsula there is an increase of the bias during June-August (Figure 20a) and in SON over the Arabian peninsula (Figure 20b). This is consistent with the already existing over-estimation of the diurnal range in the control simulation, which is likely associated with limitations of the LSA-SAF product in those regions due to Aerosols (West Sahara) and directional effects (Arabian Península). Both changes in land cover, LAI and vertical resolution in the soil show a positive impact on the LST diurnal range. Further research is needed to understand if these two are additive and to assess the benefits (and limitations) of these two changes in the representation of the LST diurnal cycle.

The impact on soil moisture is evaluated with the ESA-CCI soil moisture product (Gruber et al. [83], Dorigo et al. [84] and Gruber et al. [85]). Figure 21 shows the difference of the surface soil moisture correlation with regards to the ESA-CCI surface soil moisture product between simulations using the 10-layer configuration and the current 4-layer configuration for JJA 2018. It is shown that during dry season, the surface soil moisture of the 10-layer configuration correlates better with the ESA-CCI soil moisture as the satellite observation depth matches better the surface layer depth of the 10-layer configuration.

With the coupling with the river discharge model (Section 4.7), and the possibility brought with the hydrological verification, this development would benefit from additional constraints for uncertain or non-observable model parameters related to soil hydrology which would facilitate its operational implementation. The estimation of these uncertain model parameters would even be facilitated and automated once the multiscale parameter regionalization system (MPR [52]) is fully coupled with ECLand. This development would be then viable for operational implementation in the upcoming cycles.

4.6. Evaluation of the Updated Surface Water

Modelled lake surface temperatures were verified against in situ measurements recorded by the Data and Information Centre of the Finnish Environment Institute (SYKE) over 27 lakes for 2010–2014. In situ lake water surface temperatures are measured at 8 a.m. local time near the shore during the ice-free season 20 cm below the water surface [53]. The model values from offline surface experiments forced with ERA5 were sampled to correspond to the measured values (i.e., during ice-free season at 5:00 or 6:00 UTC depending on daylight saving time).

An overall error reduction is seen in Figure 22. With a more realistic lake depth, the mean absolute error from 2010 to 2014 becomes smaller for nearly all lakes. Averaged over all 27 lakes, a 13.4% error reduction is achieved. In some cases (Kilpisjarvi, Pesiojarvi, Inarijarvi, Rehja-Nuas, Paijanne, Kuivajarvi, Pyhajarvi and Paajarvi2), the modelling errors for lake water temperature did not change much since their depth remained close to the control experiment depth. For Lappajarvi, the temperature scores showed a deterioration because its updated depth was 2.5 times bigger than the observed values due to the depth horizontal interpolation resolution.

The updated land/water maps are planned for implementation in the ECMWF IFS operational system, as it will reduce uncertainties related to such ancillaries. In this case the evaluation is performed based on data assimilation experiments with subsequent daily forecasts. The near surface temperatures are shown in Figure 23. The simulations were performed for 2020 winter and 2019 summer seasons. RMS errors of daily 24 and 36 h forecasts are computed with respect to the own analysis at the corresponding verification time. Figure 23 shows the 2 m temperature RMSE difference, normalised by rmse, between the forecast simulations based on the new land/water maps and simulations based on the current operational maps. An overall neutral signal is seen globally for both day and night ranges (t + 36 and t + 48) during the northern hemisphere Winter season (figure not shown). Summer is characterised with a neutral to positive impact (Figure 23). Locally, a slightly negative impact is seen over the Canadian Northern coast and a positive signal marked with a reduction of the 2 m temperature RMSE when using the updated land/water maps over Greenland. These impacts correspond to the differences of land/water points seen in Figure 7. The overall evaluation shows that the IFS would benefit from a more accurate surface water representation and it is therefore considered for operational implementation in the next cycle.

4.7. Evaluation of Coupled River Discharge

The evaluation of river discharge in different configurations is presented in Figure 24. The main goal of these results is to illustrate the different capabilities of the system for various applications. The evaluation displays the distribution of the temporal correlation and percent bias of daily river discharge for several stations on the Mackenzie (panels a,b) and Danube (panel c,d) basins, with the number of stations listed on the top of each panel. The observations are taken from the Global Runoff Data Centre (GRDC). All simulations were carried out at the $1 / 4^{\circ}$ resolution for CaMa-Flood and TL639 for ECLand (same resolution as ERA5, about 32 km) extending since 1979 to 2018, and the evaluation is performed between 1982 and 2018 considering the available station data with at least 5 years of daily data. The boxplots in each panel show the distribution of the correlation and percent bias for six configurations. “coup1” is the default 1-way coupling with ECLand driven by ERA5 near-surface meteorology. “coup2” is a 2-way coupling simulation allowing water extraction from floodplains and update of inland water fraction in ECLand following CaMa-Flood flooded area extent. “iner” and “kine” are two stand-alone CaMa-Flood simulations driven by the runoff of “coup1”, differing in the flow velocity calculations in CaMa-Flood, using the local inertia in “iner” (the default configuration) and the kinematic wave in “kine”. Finally, “era5” and “era5l” are also stand-alone CaMa-Flood simulations driven by the runoff of the reanalysis ERA5 and ERA5-Land, respectively.

The results in Figure 24 show a neutral impact on river discharge simulation of the 2-way coupling, when compared with 1-way coupling on the selected basins. This is expected as both basins do not have large flooded areas. Results on some stations in the Amazon basin show some impact but reduced (not shown). However, the 2-way coupling configuration is still experimental and further developments are necessary for a consistent 2-way coupling, in particular on the handling of grid-scale versus sub-grid scale inland water in ECLand. The CaMa-Flood stand-alone configuration iner, shows identical results to coup1, as expected. The kinematic solution for flow velocities (kine) has some impact on the temporal correlations, with mixed results in the two analyzed basins. The local inertia solution is more accurate for large upstream areas where small water level slopes dominate the river flow, and is of particular relevance in representing surface water elevation. In smaller upstream regions with higher terrain slopes both local inertial and kinematic wave solutions are similar, with a reduced computational cost for the kinematic wave solution. Therefore, for particular applications were computational cost is relevant (e.g., ensemble forecasting) this could be considered. Finally, the flexibility of driving the hydrodynamic component with different runoff inputs is illustrated with the results of using ERA5 and ERA5-Land (era5 and era5l simulations in Figure 24). The era5l simulations are very close to coup1 (and iner) as the main difference is the spatial resolution of the runoff fields produced by ECLand (9 km in the case of era5l and 35 km in coup1/iner). Larger differences appear in era5, in particular for the percent bias (Figure 24b,d). These are associated with the soil moisture and snow mass increments of the land data assimilation in ERA5 that changes the land surface water balance [56]. As previously mentioned, the full operational integration of CaMa-Flood in the IFS will be possible in future cycles with additional developments for distributed memory parallelism in CaMa-Flood and for coupling of river discharge with the ocean model.

4.8. Evaluation of Urban Representation

The evaluation of the single-layer urban scheme has been performed using both the IFS single-column model and ECLand. Urban parameters were first optimised with the IFS single-column model using a Gauss-Newton non-linear least square algorithm and observations of 2 m temperature from 8 European surface sites. The surface albedo, emissivity, volumetric heat content and thermal conductivity were simultaneously optimised for average urban roofs, roads and walls. Monthly comparisons with hourly output for independent (January) observations show model improvement when using the urban scheme at all sites except Warsaw (Figure 25). The Warsaw site is at an airport located outside of the city, which does not provide a realistic representation of an urban environment. When using the urban scheme, the largest temperature changes occur in summer at night with smaller changes in the daytime (figure not shown). These changes typically agree well with observation based studies, which also find urban heat island enhancements are larger at night than during daytime [93,94].

High-resolution global ECLand simulations (∼1 km horizontal) show expected skin and 2 m temperature enhancements over both large and small urban conurbations (Figure 26). At coarser resolution the signal from smaller conurbations is lost (∼9 km) and a small urban signal averaged over an unrealistically large area is simulated at even coarser resolutions (∼30 km). It is therefore recommended that a reasonable land surface representation using the single-layer urban canopy model is only derived at resolutions of 9 km or higher. This development shows promising results to be considered as operationally viable for future cycles.

5. Summary and Perspectives

The ECLand land-surface modelling system aims to expedite progress in Earth system modelling and greatly benefits from international collaboration by co-developing and facilitating modular extensions. ECLand is based on the Carbon-Hydrology Tiled Scheme for Surface Exchanges over Land in its current operational version (CY47R1) and represents an integral part of the IFS for ECMWF’s current and future global weather and environmental (services) applications.

This paper highlights examples of recent developments in different components of the land surface model, namely: (i) vegetation, (ii) snow, (iii) soil, (iv) open water/lake, (v) river/inundation, (vi) urban areas, which are illustrated and evaluated against benchmark observational datasets. The vegetation results show potential NWP improvements with updated LU/LC maps as well as when including a dynamic vegetation representation which also benefits the accuracy of the land biogenic fluxes. The open water and lake representation has improved in accuracy thanks the most up-to-date satellite based products and with an overall neutral impact for NWP. The coupling with the river/inundation component showed potential to further improve the water representation and its variability however for an integration with the IFS, CaMa-Flood still need additional technical developments for distributed memory parallelism. The updated soil developments proved that a refined and extended layers discretisation has potential to improve the coupling with the atmosphere and the match with the soil moisture satellite observation. Additional investigation using the multiscale parameter regionalization system (MPR [52]) to improve and better optimise the soil parameters will boost future operational implementation. on the other hand, the new snow representation is confirmed as part of the next operational cycle operational as the extended evaluation proved that it improves the snow depth and reduces the 2 m temperature bias. The urban parametrization showed promising results in improving the 2 m temperature over city areas and will be further refined with global forecast evaluation. Future developments on the seamlessness of the land surface system target the improvement of processes affecting also the seasonal range and reanalyses, for example, through consideration of the vegetation status inter-annual variability as also depicted in [26].

Perspectives for NWP are focusing on improvements of the water cycle, with particular attention to evaporation processes, transition seasons, hydrological forecast skill and with a focused attention to a better exploitation of emerging very-high resolution observational datasets (e.g., SENTINEL), skin temperature over land with consideration to its compensating error-free interpretation. On the snow and soil perspectives, the multi-layer approach available in ECLand for snow and soil moisture opens avenues for enhanced forward modelling of satellite radiances which will enable developments of coupled assimilation of surface sensitive observations. Preliminary investigations were conducted by Hirahara et al. [95] showing that low frequency microwave emissivity is better represented when a multilayer snowpack model is used.

In the near future, the very high resolution modelling capabilities available with ECLand are planned to be extended to land data assimilation, using Sentinel satellites for soil moisture and snow analysis, as well as skin temperature and vegetation variables, for land reanalysis and NWP applications. Land-atmosphere coupling development priorities have been defined [26] along with ways to verify model improvements using super-sites, as well as the standard observations used in ECMWF verification activities [96]. Benchmarking activity is essential for model improvement and would also include comparison against other land surface models such as in Terzago et al. [97] or the ones under the iLAMB [98] or PLUMBER phase2 [21,22] frameworks. In addition, coupling with the river discharge model will add further constraints on uncertain model parameters and would enable investigation of coupled assimilation of inland water level information in ECLand for the benefit of river discharge and land surface processes modelling. Preliminary work using the uncoupled river/IFS system showed that SMOS data assimilation has an impact on river discharge forecasts [99]. The nature of physically based models is by default attached to uncertainties coming from non-fully adequate equations for different reasons such as scale representativeness or even non-observable parameters. Ongoing collaborative efforts in ECLand are also focused on the refactoring of the land surface model to externalize process-sensitive parameters and attach the MPR system [52] to the offline system in order to enable more flexible parameter estimation and optimization. AI/Machine learning techniques could also be considered in this context.

Finally, enhanced vegetation and photosynthesis modelling, as well as very-high resolution modelling including anthropogenic elements such as urban, CO $_{2}$ and CH $_{4}$ , and prognostic fire modelling are important elements of ECLand extensions to support current and future Copernicus services.

Future developments also focus on preserving the seamlessness of the land surface system through the improvement of the processes affecting the extended range and reanalysis, such as the consideration of the vegetation status and its inter-annual variability for upcoming C3S seasonal and reanalysis systems. This effort would be consolidated by the foreseen improved estimation and characterisation of the carbon cycle components especially the biogenic, contributing to the CAMS services.

Supplementary Materials

The Supplementary Materials are available online at https://www.mdpi.com/article/10.3390/atmos12060723/s1.

Author Contributions

Conceptualization S.B. and G.B.; methodology, formal analysis and writing original draft preparation, S.B., G.B., N.W., G.A., E.D., J.M., M.C., A.A.-P., A.B., I.H., writing review and editing, all authors. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We would like to extend our thanks to the pioneers and the different contributors to the ECMWF land surface modeling efforts: Pedro Viterbo, Bart van den Hurk, Andrea Manrique-Sunen, Annika Nordbo, Clement Albergel, Lionel Jarlan, Sebastien Lafont, Jean Christophe Calvet, Isabel Trigo, Cor Jacobs, Andrea Alessandri, Anne Verhoef, Paul Dirmeyer, Rui Salgado, Victor Stepanenko and Dmitrii Mironov for their valuable development contributions. We thank Andy Brown, Florian Pappenberger, Richard Engelen for their support and fruitful comments on the manuscript. We also thank Anabel Bowen for her help improving the figures. We would also like to thank the anonymous reviewers who made valuable suggestions that led to improve the paper. This work used eddy covariance data acquired by the FLUXNET community and in particular the Hu-mat site data provided through the ImagineS FP7 project. The WOFOST data was also provided through the ImagineS FP7 project by the European Joint research Centre, We thank in particular the PIs and researchers that decided to share their data freely within the scientific community.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

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

Figures and Table

View Image - Figure 1. Simplified representation of the main ECLand system components and its prospective evolution (inner light green box) as part of the full Integrated Forecasting System (outer yellow box) where the current ECLand main components are within the light grey box.

Figure 1. Simplified representation of the main ECLand system components and its prospective evolution (inner light green box) as part of the full Integrated Forecasting System (outer yellow box) where the current ECLand main components are within the light grey box.

View Image - Figure 2. Vegetation cover difference: ESA-CCI - GLCCv1.2 based maps: (a) for the high vegetation types (forest and trees); (b) for the low vegetation types (crops, grass and herbaceous types).

Figure 2. Vegetation cover difference: ESA-CCI - GLCCv1.2 based maps: (a) for the high vegetation types (forest and trees); (b) for the low vegetation types (crops, grass and herbaceous types).

View Image - Figure 3. LAI differences for the month of July: Conservative disaggregation—Current: (a) for the high vegetation types (forest and trees); (b) for the low vegetation types (crops, grass and herbaceous types).

Figure 3. LAI differences for the month of July: Conservative disaggregation—Current: (a) for the high vegetation types (forest and trees); (b) for the low vegetation types (crops, grass and herbaceous types).

$View Image - Figure 4. The effect of seasonal varying vegetation cover. (a) Bareground fraction difference between the seasonally varying and the static vegetation map for July. (b) The seasonal amplitude of low vegetation cover as difference between July and January.$

Figure 4. The effect of seasonal varying vegetation cover. (a) Bareground fraction difference between the seasonally varying and the static vegetation map for July. (b) The seasonal amplitude of low vegetation cover as difference between July and January.

View Image - Figure 5. Schematic of the multi-layer snow scheme (a) for a 1-metre thick snowpack with its vertical discretization (d1,⋯,d5, in m) for flat terrains. (b) Idealized time-series of snow depth accumulation and ablation (continuous line) over flat terrain, with the vertical discretization used in the multi-layer scheme. The list of symbols used in (a,b) is also reported.

Figure 5. Schematic of the multi-layer snow scheme (a) for a 1-metre thick snowpack with its vertical discretization (d1,⋯,d5, in m) for flat terrains. (b) Idealized time-series of snow depth accumulation and ablation (continuous line) over flat terrain, with the vertical discretization used in the multi-layer scheme. The list of symbols used in (a,b) is also reported.

View Image - Figure 6. Updated soil depth and vertical discretization. (Left) current operational configuration with 2.89 m depth and 4 layers (7 cm, 21 cm, 72 cm, 189 cm). (Right) the updated configuration with 8 m depth and 10 layers (1 cm, 2 cm, 4 cm, 9 cm, 12 cm, 30 cm, 42 cm, 100 cm, 200 cm, 400 cm).

Figure 6. Updated soil depth and vertical discretization. (Left) current operational configuration with 2.89 m depth and 4 layers (7 cm, 21 cm, 72 cm, 189 cm). (Right) the updated configuration with 8 m depth and 10 layers (1 cm, 2 cm, 4 cm, 9 cm, 12 cm, 30 cm, 42 cm, 100 cm, 200 cm, 400 cm).

Figure 7. Lake cover difference: new GSWE based map—operational map.

View Image - Figure 8. Schematic of the single-layer urban model used in ECLand including key parameters (α: albedo, ε: emissivity, λ: thermal conductivity, C: volumetric heat capacity), variables (T: temperature, Q: moisture content), and fluxes of heat and moisture (ΦH,W).

Figure 8. Schematic of the single-layer urban model used in ECLand including key parameters (α: albedo, ε: emissivity, λ: thermal conductivity, C: volumetric heat capacity), variables (T: temperature, Q: moisture content), and fluxes of heat and moisture (ΦH,W).

View Image - Figure 9. Mean diurnal range (drange) bias of Land Surface Temperature (LST) in the control simulation (CTR) during June-August (JJA, (a)) and September-November (SON, (b)). Difference of the mean absolute bias in the simulation with new ESA-CCI LU/LA map and LAI disagragation (CLAI) in respect to the control simulation for JJA (c) and SON (d). Note that the color bar for differences in (c,d) has a range that is only half of the range in (a,b). The statistics were computed for the period 2004 to 2015 considering only clear-sky conditions.

Figure 9. Mean diurnal range (drange) bias of Land Surface Temperature (LST) in the control simulation (CTR) during June-August (JJA, (a)) and September-November (SON, (b)). Difference of the mean absolute bias in the simulation with new ESA-CCI LU/LA map and LAI disagragation (CLAI) in respect to the control simulation for JJA (c) and SON (d). Note that the color bar for differences in (c,d) has a range that is only half of the range in (a,b). The statistics were computed for the period 2004 to 2015 considering only clear-sky conditions.

View Image - Figure 10. Annual mean evapotranspiration relative bias (%) with respect to GLEAM (a) and FluxCom (c); and the absolute relative bias difference (%) between the ESA-CCI + LAI (CLAI) experiment and the control experiment (CTR), where GLEAM is used for verification in (b) and FluxCom in (d). Note the different range of the color-bars between (a,c) (−100 to 100) and (b,d) (−50 to 50).

Figure 10. Annual mean evapotranspiration relative bias (%) with respect to GLEAM (a) and FluxCom (c); and the absolute relative bias difference (%) between the ESA-CCI + LAI (CLAI) experiment and the control experiment (CTR), where GLEAM is used for verification in (b) and FluxCom in (d). Note the different range of the color-bars between (a,c) (−100 to 100) and (b,d) (−50 to 50).

View Image - Figure 11. A 2 m temperature normalised RMSE difference between 36-h forecasts using ESA-CCI LU/LC with the new LAI disaggregation operator (CLAI) and a control simulation using the operational configuration (CTR) for JF (a), MAM (b), JJA (c) and SON (d).

Figure 11. A 2 m temperature normalised RMSE difference between 36-h forecasts using ESA-CCI LU/LC with the new LAI disaggregation operator (CLAI) and a control simulation using the operational configuration (CTR) for JF (a), MAM (b), JJA (c) and SON (d).

View Image - Figure 12. 2 m temperature normalised RMSE difference between forecasts using the seasonal vegetation cover (FSCV) and the control simulations with static vegetation cover (CTR) for (a) DJF, (b) MAM, (c) JJA and (d) SON.

Figure 12. 2 m temperature normalised RMSE difference between forecasts using the seasonal vegetation cover (FSCV) and the control simulations with static vegetation cover (CTR) for (a) DJF, (b) MAM, (c) JJA and (d) SON.

Figure 13. Above ground biomass-based Anomaly Index (AIAGB) for November 2010 in % of the 1999–2013 mean.

View Image - Figure 14. In-Situ evaluation of the vegetation related variables when using the growth model. (a) Observed LAI in m2/m2 for 2004 at Fluxnet Hu-mat station (green with red dots), LAI climatology (purple), and prognostic LAI (black); (b) Above-ground biomass in kg/m2 from WOFOST (red), ECLand with climatological LAI (green) and ECLand with prognostic LAI (black).

Figure 14. In-Situ evaluation of the vegetation related variables when using the growth model. (a) Observed LAI in m2/m2 for 2004 at Fluxnet Hu-mat station (green with red dots), LAI climatology (purple), and prognostic LAI (black); (b) Above-ground biomass in kg/m2 from WOFOST (red), ECLand with climatological LAI (green) and ECLand with prognostic LAI (black).

View Image - Figure 15. Taylor diagram of Snow Water Equivalent (SWE, left) and snow depth (SD, right) for offline simulations using the single-layer (CTL SL, red) and multi-layer (ML5, blue) snow schemes over the ESM-SnowMIP sites; see Section 4.2 for details on the ESM-SnowMIP sites.

Figure 15. Taylor diagram of Snow Water Equivalent (SWE, left) and snow depth (SD, right) for offline simulations using the single-layer (CTL SL, red) and multi-layer (ML5, blue) snow schemes over the ESM-SnowMIP sites; see Section 4.2 for details on the ESM-SnowMIP sites.

View Image - Figure 16. Time series of snow water equivalent (SWE), snow Depth, and snow Density for Swamp Angel (swa) from observations (Obs, black dots) and the offline experiments using the single-layer (SL, red lines), multilayer (ML5, blue lines), and multilayer without the snow depth-dependent snow layer thickness described in Section 3.2 (ML-noOro, green lines). The latter configuration corresponds to a simulation with the same vertical discretization of the snowpack as used by Arduini et al. [41].

Figure 16. Time series of snow water equivalent (SWE), snow Depth, and snow Density for Swamp Angel (swa) from observations (Obs, black dots) and the offline experiments using the single-layer (SL, red lines), multilayer (ML5, blue lines), and multilayer without the snow depth-dependent snow layer thickness described in Section 3.2 (ML-noOro, green lines). The latter configuration corresponds to a simulation with the same vertical discretization of the snowpack as used by Arduini et al. [41].

View Image - Figure 17. (Left) RMSE difference of snow depth between global offline simulation using the multi-layer (ML5) and single-layer (SL) snow schemes for December to May period from 2014 to 2018. Blue colours indicate reduced RMSE in ML5. (Right) boxplot of snow depth statistics (bias and RMSE) using stations considered in left panel.

Figure 17. (Left) RMSE difference of snow depth between global offline simulation using the multi-layer (ML5) and single-layer (SL) snow schemes for December to May period from 2014 to 2018. Blue colours indicate reduced RMSE in ML5. (Right) boxplot of snow depth statistics (bias and RMSE) using stations considered in left panel.

View Image - Figure 18. Temperature verification at screen level (T2m). (a) Absolute mean error difference (amef) at a lead time of 2 days between coupled forecasts with the multi-layer (snowML) and single-layer (snowSL) snow schemes; (b) Same as (a), but showing root-mean-square-error difference (rmsef).

Figure 18. Temperature verification at screen level (T2m). (a) Absolute mean error difference (amef) at a lead time of 2 days between coupled forecasts with the multi-layer (snowML) and single-layer (snowSL) snow schemes; (b) Same as (a), but showing root-mean-square-error difference (rmsef).

View Image - Figure 19. Difference in correlation with FluxCom energy fluxes between the 10-layers soil experiment (OSML10_TW) and the standard 4-layers experiment (a) Sensible heat flux; (b) Latent heat flux.

Figure 19. Difference in correlation with FluxCom energy fluxes between the 10-layers soil experiment (OSML10_TW) and the standard 4-layers experiment (a) Sensible heat flux; (b) Latent heat flux.

View Image - Figure 20. Difference in absolute mean diurnal range (drange) bias of Land Surface Temperature (LST) between the simulation with OSML10_TW (10L) and the control (CTR) simulation for JJA (a) and SON (b). The statistics were computed for the period 2004 to 2015 considering only clear-sky conditions. The biases of the control simulation are shown Figure 9.

Figure 20. Difference in absolute mean diurnal range (drange) bias of Land Surface Temperature (LST) between the simulation with OSML10_TW (10L) and the control (CTR) simulation for JJA (a) and SON (b). The statistics were computed for the period 2004 to 2015 considering only clear-sky conditions. The biases of the control simulation are shown Figure 9.

View Image - Figure 21. Difference of the surface soil moisture correlation with regards to the ESA-CCI surface soil moisture product between simulations using the 10-layer configuration and the current 4-layer configuration for JJA.

Figure 21. Difference of the surface soil moisture correlation with regards to the ESA-CCI surface soil moisture product between simulations using the 10-layer configuration and the current 4-layer configuration for JJA.

View Image - Figure 22. Lake temperature mean absolute error for 27 in-situ Finnish lakes over the 2010–2014 period. Results for the current operational system are in red and results with the updated lake depths are in green.

Figure 22. Lake temperature mean absolute error for 27 in-situ Finnish lakes over the 2010–2014 period. Results for the current operational system are in red and results with the updated lake depths are in green.

View Image - Figure 23. A 2 m temperature normalised root mean square error (RMSE) difference for northern hemisphere Summer between 4dvar simulations with the updated surface water map and the current operational version. The error differences refer to 36-h (a) and 48-h (b) forecasts initialized at 0UTC.

Figure 23. A 2 m temperature normalised root mean square error (RMSE) difference for northern hemisphere Summer between 4dvar simulations with the updated surface water map and the current operational version. The error differences refer to 36-h (a) and 48-h (b) forecasts initialized at 0UTC.

View Image - Figure 24. River discharge evaluation distribution of several stations in the Mackenzie (a,b) and Danube (c,d) basins. The evaluation is performed on daily data, displaying the distribution of correlation (a,c) and percent bias (b,d) of several model configurations (vertical axis): coup1: 1-way coupling between ECLand and CaMa-Flood; coup2: 2-way coupling between ECLand and CaMa-Flood stand-alone simulations with the flow velocity calculations using the local inertia solution (iner) and kinematic wave (kine), both driven by runoff from coup1. era5 and era5l are also stand-alone CaMa-Flood simulations driven by the runoff of ERA5 and ERA5-Land reanalysis, respectively. The boxplots whiskers extend from percentiles 5 to 95, with the box representing percentiles 25 and 75 and the median in red. The values displayed on the left, close to the vertical axis, identify the median value for each distribution.

Figure 24. River discharge evaluation distribution of several stations in the Mackenzie (a,b) and Danube (c,d) basins. The evaluation is performed on daily data, displaying the distribution of correlation (a,c) and percent bias (b,d) of several model configurations (vertical axis): coup1: 1-way coupling between ECLand and CaMa-Flood; coup2: 2-way coupling between ECLand and CaMa-Flood stand-alone simulations with the flow velocity calculations using the local inertia solution (iner) and kinematic wave (kine), both driven by runoff from coup1. era5 and era5l are also stand-alone CaMa-Flood simulations driven by the runoff of ERA5 and ERA5-Land reanalysis, respectively. The boxplots whiskers extend from percentiles 5 to 95, with the box representing percentiles 25 and 75 and the median in red. The values displayed on the left, close to the vertical axis, identify the median value for each distribution.

View Image - Figure 25. Hourly observed 2 m temperature taken from 8 urban sites for January 2012 (black circles); also shown are single-column model results using the control model (SCM_43R3: red) and the optimised urban scheme (SCM-URB: blue). Numbers indicate RMSE values when compared to observations (adapted from [69]).

Figure 25. Hourly observed 2 m temperature taken from 8 urban sites for January 2012 (black circles); also shown are single-column model results using the control model (SCM_43R3: red) and the optimised urban scheme (SCM-URB: blue). Numbers indicate RMSE values when compared to observations (adapted from [69]).

View Image - Figure 26. The difference between urban and control monthly mean 2 m temperature over central UK for 00:00 UTC January 2019 at 30 km (top), 9 km (middle) and 1 km horizontal resolution (bottom) taken from ECLand simulations. Conurbations larger than 1000 km2 (solid), 500 km2 (dashed) and 100 km2 (dotted) are denoted with boxes (adapted from [69]).

Figure 26. The difference between urban and control monthly mean 2 m temperature over central UK for 00:00 UTC January 2019 at 30 km (top), 9 km (middle) and 1 km horizontal resolution (bottom) taken from ECLand simulations. Conurbations larger than 1000 km2 (solid), 500 km2 (dashed) and 100 km2 (dotted) are denoted with boxes (adapted from [69]).

Table 1

Experimental configurations for the evaluation of model developments.

Model Development	Experiment Type	Name	Comment
Surface water representation	Offline Surface simulation	OSW	Updated lake parameters
Surface water representation	Coupled Forecast + Data assimilation	DASW	Updated Glacier, Land water mask and lake
Updated LU/LC and LAI	Offline Surface simulation	CLAI	ESA-CCI land cover and LAI conservative disaggregation
Updated LU/LC and LAI	Coupled Forecast simulation	FCLAI	ESA-CCI land cover and LAI conservative disaggregation
Seasonal vegetation cover	Offline Surface simulation	OSCV	LAI based Seasonal vegetation cover
Seasonal vegetation cover	Coupled Forecast simulation	FSCV	LAI based Seasonal vegetation cover
Prognostic LAI	Offline Surface simulation	OPLAI	Prognostic LAI based on daily CO $_{2}$ balance
New snow multilayers	Offline Surface simulation	OSnML5	5 Layers snow scheme
New snow multilayers	Coupled Forecast	FSnML5	5 layers scheme
Increased soil depth and discretization	Offline Surface simulation	OSML10_TW	10 Layers soil 8 m for temperature and water
Increased soil depth and discretization	Offline Surface simulation	OSML10_T	10 Layers soil 8 m for temperature and 3 m for water
Coupling with hydrodynamic	Offline Surface simulation coupled with routing model	-	details in Section 4.7
Urban parametrization	Offline Surface simulation	OSMU	offline 1D and 2D urban parametrization
Urban parametrization	Single Column simulation	SCMU	coupled 1D forecast with urban parametrization

Word count: 13991

Show less

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Abstract

Translate

The land-surface developments of the European Centre for Medium-range Weather Forecasts (ECMWF) are based on the Carbon-Hydrology Tiled Scheme for Surface Exchanges over Land (CHTESSEL) and form an integral part of the Integrated Forecasting System (IFS), supporting a wide range of global weather, climate and environmental applications. In order to structure, coordinate and focus future developments and benefit from international collaboration in new areas, a flexible system named ECLand, which would facilitate modular extensions to support numerical weather prediction (NWP) and society-relevant operational services, for example, Copernicus, is presented. This paper introduces recent examples of novel ECLand developments on (i) vegetation; (ii) snow; (iii) soil; (iv) open water/lake; (v) river/inundation; and (vi) urban areas. The developments are evaluated separately with long-range, atmosphere-forced surface offline simulations and coupled land-atmosphere-ocean experiments. This illustrates the benchmark criteria for assessing both process fidelity with regards to land surface fluxes and reservoirs of the water-energy-carbon exchange on the one hand, and on the other hand the requirements of ECMWF’s NWP, climate and atmospheric composition monitoring services using an Earth system assimilation and prediction framework.

Details

Title

ECLand: The ECMWF Land Surface Modelling System

Author

Boussetta, Souhail¹; Balsamo, Gianpaolo¹

; Arduini, Gabriele¹; Dutra, Emanuel²

; McNorton, Joe¹; Choulga, Margarita¹; Agustí-Panareda, Anna¹; Beljaars, Anton¹; Wedi, Nils¹; Munõz-Sabater, Joaquín¹

; de Rosnay, Patricia¹

; Sandu, Irina¹; Hadade, Ioan¹; Carver, Glenn¹

; Mazzetti, Cinzia¹; Prudhomme, Christel¹; Yamazaki, Dai³; Ervin Zsoter¹

¹ European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG2 9AX, UK; [email protected] (G.A.); [email protected] (J.M.); [email protected] (M.C.); [email protected] (A.A.-P.); [email protected] (A.B.); [email protected] (N.W.); [email protected] (J.M.-S.); [email protected] (P.d.R.); [email protected] (I.S.); [email protected] (I.H.); [email protected] (G.C.); [email protected] (C.M.); [email protected] (C.P.); [email protected] (E.Z.)
² Instituto Português do Mar e da Atmosfera, 1749-077 Lisbon, Portugal; [email protected]; Instituto Dom Luiz, IDL, Faculty of Sciences, University of Lisbon, 1749-016 Lisbon, Portugal
³ Institute of Industrial Science, The University of Tokyo, Tokyo 153-8505, Japan; [email protected]

First page

723

Publication year

2021

Publication date

2021

Publisher

MDPI AG

e-ISSN

20734433

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.3390/atmos12060723

ProQuest document ID

2544577401

ECLand: The ECMWF Land Surface Modelling System

Jump to:

Full text

Abstract

Details

Suggested sources