ECLand: The ECMWF Land surface modelling system

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


Introduction
values of k range between 0.4 and 0.8, here we consider a value of 0.6 which was also adopted by 224 Nogueira et al. [37]. In the latter study, when the clumping was also applied to the high vegetation, 225 coupled simulations indicated a strong impact over the Northern Hemisphere circulation which require 226 further investigations, therefore the implementation in this study is so far limited to low vegetation. a) Bareground fraction difference between seasonal varying and static vegetation (July).
b) Low vegetation cover difference between July and January where N l is the leaf Nitrogen concentration, and e, f are respectively the slope and intercept nitrogen where r B is the relaxation coefficient. For the sake of completeness, the growth model is further detailed 240 in the supplementary materials.

252
The heat budget of the snowpack is solved implicitly for the active snow layers, and using an 253 explicit flux coupling at the snow-atmosphere interface [42]. To ensure numerical stability over forested (first-guess of the snow temperature). Movement of liquid water between layers (meltwater or rainfall) 261 is taken into account using a bucket-type formulation as in [16]. The latent heat released or absorbed 262 during phase changes in a snow layer is used to vary the temperature of the respective snow layer.

263
The number of active snow layers and their thicknesses is computed diagnostically at the 264 beginning of each time-step before the prognostic snow fields are updated. resolution at the interfaces to be maintained with the atmosphere above and the soil underneath. An 274 example of the vertical discretization of a 1.0 m-thick snowpack is shown in Fig. 5c.

275
Over complex terrains, for D sn < 0.25 m the same vertical discretization as described for a flat terrain is used. For a snowpack with D sn ≥ 0.25 m, the minimum (D min,i ) and maximum (D max,i ) thicknesses of each layer (i.e. the minimum and maximum allowed thickness for an active snow layer) varies with snow depth as follows:   In addition to the structural aspects described above, ML differs from SL also in the   With higher spatial resolution and more accurate grid parameters characterisation, in addition to 294 broader applications and the need to improve the water budget in the model and coupled hydrological 295 components, the importance of the soil depth and vertical discretisation within land surface models has 296 increased. It affects the way land surfaces store and regulate water, energy and also carbon fluxes. The 297 amount of water in the soil and its vertical distribution in the column are important for the regulation 298 of heat and water vapour fluxes towards the atmosphere spanning a range of time scales from minutes 299 to months in the coupled land-atmosphere system. This is further modulated by the vertical roots 300 distribution, and soil moisture stress function, which control evapotranspiration under soil moisture 301 stress conditions (next section). Currently in the ECMWF land Surface Scheme the soil column is 302 represented by a fixed 4-layer configuration with a total of 2.89m depth. In the ECLand system we 303 introduce new configurations with increased soil depth (up to 8m) and higher vertical discretisation  An optional formulation for the distribution of roots density has been also developed in ECLand.

308
It is based on a uniform distribution, which is appealing due to its simplicity and straightforward

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

325
The upgraded inland water mask is based on the upgraded fractional land-water mask at 1 km 326 resolution. First, fractional land-water mask is translated into "land" and "water" binary mask, then by 327 using an innovative "waterbody" separation algorithm [53] "water" is separated into ocean and inland  The key characteristic of CaMa-Flood is the explicit representation of flood stage, representing 366 water level and inundation extent, which complements the conventional use of river discharge [62].  CaMa-Flood source code is integrated in the ECLand system allowing stand-alone CaMa-Flood 371 simulations or simulations with 1-way and 2-way coupling between ECLand and CaMa-Flood.

372
The coupling is performed sequentially and the code is integrated as a single executable (without 373 external couplers) following a similar approach as the coupling with the ocean model to IFS [68]. In 374 1-way coupling, the default configuration, surface and sub-surface runoff are sent to CaMa-Flood.

375
In stand-alone mode CaMa-Flood reads from input files the runoff data, reproducing the 1-way 376 coupling results. This configuration has a lower computational cost, allowing for testing of the 377 hydrodynamic model as well as driving the model with other runoff data sources like reanalysis.

378
The 2-way coupling between ECLand and CaMa-Flood was also experimentally developed allowing 379 initial testing of the interaction between inland water in ECLand and flooded areas in CaMa-Flood. 380 Inland water evaporation (in ECLand) can be extracted from flooded areas in CaMa-Flood, and the 381 flooded area extent can update the inland water fraction in ECLand. CaMa-Flood integration in the 382 ECLand system is only available in offline mode supporting shared memory parallelism via OpenMP.

383
A full integration of CaMa-Flood in the IFS is envisaged in upcoming cycles to support coupled 384 atmosphere-ocean-land-rivers simulations aiming at a full representation of the water cycle. This will 385 still require developments both at the level of distributed memory parallelism in CaMa-Flood and in 386 the coupling of river discharge with the ocean model.

387
The river network description used in CaMa-Flood is currently available in several global regular 388 latitude/longitude resolutions: 1/4 • , 1/10 • , 1/12 • , 1/20 • , and 1/60 • . However, each unit-catchment 389 in CaMa-Flood has an irregular shape and ECLand can be setup on any regular or the IFS model's 390 Gaussian reduced grid. To accommodate the different spatial grids between these two components, 391 the fluxes are interpolated using pre-computed weights that account for the area of interception of 392 each grid element. For the 2-way coupling, the inverse of the weights is used to guarantee mass 393 conservation. These weights are computed as part of the pre-processing of the system, which also 394 allows the definition of regional domains, which are relevant for testing and for applications that do where C sn , C U dry and C wet represent the snow covered, dry urban and wet fraction, respectively.

419
For land surface skin, future developments will merge the urban cover product with the remaining 420 land cover so this re-weighting process is no longer required. For certain elements, e.g. hydrology, this 421 prevents a weighted average that neither fits an urban or non-urban regime. Currently, for sub-surface

Model efficiency 426
The surface model inherits the same "array-block" data structures as all other IFS physics 427 parameterisation packages and is therefore capable of benefiting from modern processors and GPUs.   The new developments described above are evaluated against the following observations or   Screen-level temperature and dewpoint can be verified using the analysis, which draws closely to 592 the surface synoptic observations. Figure 11 shows the 2m temperature normalised RMSE difference   Results indicate a reduction of the mean error in particular in the high-latitude regions, which 706 can be associated to an improved prediction of very cold temperatures over snow-covered regions.

707
The RMSE of T2m is reduced in the snow multi-layer forecasts over most of Eurasia regions, however 708 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 710 which the new scheme is more exposed [more "responsive", as described in 41]. However, the increase 711 of RMSE over complex terrain can be also related to an increase of the preexisting biases in the control 712 forecasts, which can be partly due to the relatively low horizontal resolution used in these experiments   Similarly to the reduction of the diurnal range underestimation seen in CLAI (Figure 9), the increased The impact on soil moisture is evaluated with the ESA-CCI soil moisture product ( Gruber et al.  The evaluation of river discharge in different configurations is presented in Figure 24. The main 792 goal of these results is to illustrate the different capabilities of the system for various applications.

793
The evaluation displays the distribution of the temporal correlation and percent bias of daily river configuration) and the kinematic wave in "kine". Finally, "era5" and "era5l" are also stand-alone 806 CaMa-Flood simulations driven by the runoff of the reanalysis ERA5 and ERA5-Land, respectively.

807
The results in Figure 24 show a neutral impact on river discharge simulation of the 2-way coupling,

850
The ECLand land-surface modelling system aims to expedite progress in Earth system modelling 851 and greatly benefits from international collaboration by co-developing and facilitating modular Perspectives for NWP are focusing on improvements of the water cycle, with particular attention 862 to evaporation processes, transition seasons, hydrological forecast skill and with a focused attention to 863 a better exploitation of emerging very-high resolution observational datasets (eg. SENTINEL), skin 864 temperature over land with consideration to its compensating error-free interpretation. On the snow 865 and soil perspectives, the multi-layer approach available in ECLand for snow and soil moisture opens 866 avenues for enhanced forward modelling of satellite radiances which will enable developments of 867 coupled assimilation of surface sensitive observations. Preliminary investigations were conducted by 868 Hirahara et al. [96] showing that low frequency microwave emissivity is better represented when a 869 multilayer snowpack model is used.

870
In the near future, the very high resolution modelling capabilities available with ECLand are 871 planned to be extended to land data assimilation, using Sentinel satellites for soil moisture and