The Cycle 46 Configuration of the HARMONIE-AROME Forecast Model

Abstract

The aim of this technical note is to describe the Cycle 46 reference configuration of the HARMONIE-AROME convection-permitting numerical weather prediction model. HARMONIE-AROME is one of the canonical system configurations that is developed, maintained, and validated in the ACCORD consortium, a collaboration of 26 countries in Europe and northern Africa on short-range mesoscale numerical weather prediction. This technical note describes updates to the physical parametrizations, both upper-air and surface, configuration choices such as lateral boundary conditions, model levels, horizontal resolution, model time step, and databases associated with the model, such as for physiography and aerosols. Much of the physics developments are related to improving the representation of clouds in the model, including developments in the turbulence, shallow convection, and statistical cloud scheme, as well as changes in radiation and cloud microphysics concerning cloud droplet number concentration and longwave cloud liquid optical properties. Near real-time aerosols and the ICE-T microphysics scheme, which improves the representation of supercooled liquid, and a wind farm parametrization have been added as options. Surface-wise, one of the main advances is the implementation of the lake model FLake. An outlook on upcoming developments is also included.

1. Introduction

There is a very strong history of collaboration in the area of operational numerical weather prediction (NWP) among European meteorological services. Some of this history is detailed in Bengtsson et al. [1], where the High Resolution Limited Area Model (HIRLAM) international research program and the Aire Limitée Adaptation Dynamique Développement International (ALADIN) cooperation are introduced. In 2014, HIRLAM and ALADIN agreed on the formation of a single, united consortium, which came to fruition on 27 November 2020. The new consortium, called A Consortium for COnvection-scale modeling Research and Development (ACCORD), is a collaboration between the HIRLAM, ALADIN, and the Regional Co-operation for Limited Area Modeling in Central Europe (RC-LACE) consortia. It is a unique cooperation between 26 countries in Europe and Northern Africa.

ACCORD cooperates on the development of a limited area NWP system and shares software for parts of it. Historically, there are three canonical system configurations (CSCs) within the ACCORD NWP system maintained by different groupings of countries:

• AROME (Applications of Research to Operations at Mesoscale)
• HARMONIE-AROME (HIRLAM–ALADIN Research on Mesoscale Operational NWP in Euromed)
• ALARO (Aire Limitee Adaptation/Application de la Recherche a l’Operationnel).

The HARMONIE-AROME canonical configuration is developed, maintained, and validated by the HIRLAM countries (Denmark, Estonia, Finland, Iceland, Ireland, the Netherlands, Norway, Spain, Sweden, and Lithuania), as introduced in [1]. HARMONIE–AROME is a limited-area, mesoscale, spectral, nonhydrostatic model. It features a dynamical core developed within the ALADIN consortium [2,3]. However, several adaptations and improvements were made to the model by scientists from the HIRLAM countries, mainly to the atmospheric and surface physical parametrizations. The model incorporates parametrizations for shortwave and longwave radiation, land-surface processes, cloud dynamics and microphysics, turbulence, and shallow convection. With a grid spacing of 2.5 km, the model is designed to resolve and explicitly simulate deep convection through its nonhydrostatic dynamics, eliminating the need for a deep convection parametrization.

ACCORD shares some developments (parts of the code) with the European Center for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) model. This is possible due to the modularity of the code. To ensure consistency of the common code, ACCORD follows the releases of the IFS, known as cycles. When a new IFS cycle is issued, parts of the HARMONIE-AROME code are adjusted to collocate with it, and new developments are merged in. The basic cycle numbers for HARMONIE–AROME follow the IFS cycle numbering. Local developments of HARMONIE–AROME, starting from the basic IFS release, may lead to additional local releases. Because the process of merging and harmonization is time-consuming, some IFS cycles are skipped. In this paper, we present the latest release of the HARMONIE-AROME reference system, which is Cycle 46. Previous HARMONIE-AROME releases were based on Cycle 43. The release described in [1] was related to Cycle 40, and the next release will be Cycle 49.

The HARMONIE-AROME NWP system, in addition to the forecast model, includes data assimilation (DA) to initialize the model from observations, a system for the generation of external parameter fields (such as land-use), a system to provide lateral boundary conditions from a host model, post-processing, verification, and ensemble forecasting. HARMONIE-AROME is used both for research and operational purposes. The configurations used by national meteorological services (NMSs), or their operational groupings, are usually close to the official releases, although they may deviate slightly. Currently two operational groupings use the HARMONIE-AROME NWP system. The MetCoOp grouping started in 2010 with the Norwegian Meteorological Institute and the Swedish Meteorological and Hydrological Institute collaborating on the production of operational weather forecasts. The Finnish Meteorological Institute joined in 2017, with Estonia, Latvia, and Lithuania following in 2022. Another grouping, called United Weather Centers—West (UWC-West), was formed in 2021 by the NMSs of Ireland, Denmark, Iceland and the Netherlands. The NMS of Spain, AEMET, runs HARMONIE-AROME operationally using their own computational and staffing resources. In addition, the Norwegian Meteorological Institute runs it operationally over the Arctic domain.

In this paper we give an overview of the HARMONIE-AROME forecasting model developments for the 7-year period between Cycle 40 (described in [1]) and the latest Cycle 46. The DA system and ensemble system are not covered here. As well as the reference configuration settings, some research components, which will be used in the near future, are described. The work is outlined as follows: changes in dynamics and the model configuration, as well as system aspects, are described in Section 2, changes in upper-air physics in Section 3, and changes in surface physics are detailed in Section 4. In each section, we describe briefly the current status and list changes since Cycle 40. Upcoming physics developments (upper-air and surface) are discussed in Section 5. In Section 6, the general results are discussed and conclusions are given.

2. Dynamics, Model Configuration, and System Aspects

Here, a brief overview of the dynamical core of the HARMONIE-AROME system for Cycle 46 is given. The detailed description can be found in Bengtsson et al. [1]. Information about the model configuration and system are also provided.

The nonhydrostatic core of HARMONIE-AROME is shared with the ALADIN and AROME CSCs [4,5] and was developed within the ALADIN consortium [2,3]. It uses the fully compressible Euler equations with a mass-based hybrid pressure terrain-following vertical coordinate [6] and Cartesian horizontal coordinates on several suggested map projections (such as the Lambert projection). The equations are integrated on an A grid using a two-time-level semi-Lagrangian (SL) advection scheme with semi-implicit (SI) adjustment. Double Fourier decomposition is used to compute derivatives in spectral space, while finite differencing is used in the vertical [7].

The time scheme allows for a general iterative centered implicit (ICI) approach for improved stability [3]. Typically the ALADIN and AROME consortia use a predictor–corrector method in their operational configurations, while HARMONIE-AROME uses a single-step SISL(semi-Lagrangian semi-implicit) discretization. The two approaches have a comparable cost, as the predictor–corrector allows for a slightly longer time step that off-sets the additional effort. Stability is enhanced through the use of the stable extrapolation two-time-level scheme (SETTLS) [8], which is applied to the nonlinear terms as well as in SL trajectory calculations. In addition, the implementation of the boundary relaxation method of Davies [9] for the upper boundary condition has further helped to stabilize integrations.

Linear spectral horizontal diffusion is applied to the spectral prognostic variables. In addition, the flow-dependent nonlinear semi-Lagrangian horizontal diffusion (SLHD) scheme [10] is used for hydrometeors and cloud water. Similar to that described in Seity et al. [5], quasi-monotonic interpolation is used in the semi-Lagrangian advection. Moreover, the continuous mapping about departure points (COMAD) scheme [11] is employed.

The default reference configuration grid spacing of HARMONIE-AROME is 2.5 km, with the corresponding spectral grid using a linear truncation. The vertical grid consists of 65 levels, ranging from 12 m above ground at the lowest level up to a model top of 10 hPa. The default time step is 75 s. A 16th-order diffusion operator is used to smooth the input orography. The lateral boundary conditions are generally taken from the IFS HRES forecasts of ECMWF.

As described in [1], there is an option to run the forecasting model with the spectral truncation reduced to a quadratic or cubic grid. The physical grid-spacing remains unchanged, but the formal resolution is reduced with a lower wavenumber truncation and with a reduction in computational cost of up to 20%. The quadratic grid in particular has been used successfully in operational runs of some NMSs, with a minimal impact on accuracy (e.g., Met Éireann has used it since cycle 40 and UWC-West now uses it with cycle 43).

By default, HARMONIE-AROME runs in double precision (i.e., it uses 64 bit variables for storing floating-point data). However, since Cycle 43h2 (one of the cycle 43 releases), there is a possibility to run the forecast model in single precision (with 32 bit floating-point variables), keeping all other model components in double precision. With this approach, highly precision-sensitive parts, such as data assimilation, are not impacted. This approach was adopted at ECMWF and has been used by them operationally since their Cycle 47r2 [12], providing a runtime saving of approximately 40% without negatively impacting forecast quality. Limited testing of this approach in HARMONIE-AROME Cycle 46 showed runtime savings of about 35%, with little to no impact on point verification scores apart from a small positive mean sea level pressure bias. After extensive testing, single precision forecasts are now used operationally by the UWC-West grouping with HARMONIE-AROME Cycle 43.

Within ACCORD, only the FORTRAN code of the CSCs is shared, but the infrastructure to run the model in research or operational mode is not shared. To run the HARMONIE-AROME CSC, both for research and operations, a HARMONIE-specific software infrastructure is used. The main part of this is the so-called Harmonie scripting system, written in several system scripting languages, Python, perl, etc. The scripting system organizes the data flow for all parts of the system, including the generation of external parameters, the management and quality control of the observations, variational data assimilation, the preparation of boundary conditions, running the forecast, post-processing, and the extraction of data for model verification. The model can run in ensemble mode, and this possibility is also provided by the scripting system. A work-flow manager called ecFlow, developed by ECMWF, is also part of the HARMONIE software infrastructure. A schematic of the HARMONIE-AROME workflow is shown in Figure 1.

A single-column model framework, known as Modèle Unifié, Simple Colonne (MUSC) [13], is also available within the ACCORD consortium. It is useful for testing and validating new model developments and has been used as a validation tool for some of the work described further in the following sections.

3. Upper-Air Physics

3.1. Radiation

The radiation scheme developed at ECMWF and available in the IFS starting from cycle 25R1 [14] is the default shortwave (SW) radiation scheme in HARMONIE–AROME. It contains six spectral intervals. The default longwave (LW) radiation scheme uses the Rapid Radiative Transfer Model (RRTM) of Mlawer et al. [15]; it contains 16 spectral bands. Details of both schemes can be found in the IFS [14] and the mesoscale research model Meso-NH [16] documentation. Because of computational constraints, the full radiation calculations are currently performed approximately every 15 min.

An alternative broadband radiation scheme called ACRANEB2 [17,18] is also available in HARMONIE-AROME, though it is not used by default or operationally, as explained in Bengtsson et al. [1]. Comparisons of the clear-sky SW and LW radiation fluxes at the surface with observations over Finland were performed with experiments using the default radiation scheme and ACRANEB2 [19,20]. A general agreement between the schemes was found.

Changes in the radiation scheme of HARMONIE–AROME since Cycle 40 are related mainly to the effects of clouds and aerosols. Clouds and aerosols cause the biggest current uncertainties in both SW and LW radiative fluxes in atmospheric models. In particular, the effect of aerosols on clouds is generally recognized [21] as the biggest radiative forcing uncertainty. The radiative transfer calculations use optical thickness, single scattering albedo, and asymmetry factor of cloud particles, atmospheric gases, and aerosols. In the model, information about clouds is given by prognostic specific cloud liquid water and cloud ice content (in HARMONIE-AROME, part of the precipitating in-cloud snow and graupel content is added to the cloud ice content). The atmospheric gases, which are taken into account in the radiation scheme, are (i) H₂O (water vapor), the amount of which is a prognostic model variable; (ii) CO₂, N₂O, CH₄, and O₂, the amount of which is fixed; and (iii) O₃, the amount of which is represented by a monthly climatology.

In the default version of HARMONIE-AROME cycle 46, aerosols are represented by monthly climatologies. For that, the Tegen [22] aerosol climatology is used. The use of near real-time (NRT) aerosol data, provided by the Copernicus Atmosphere Monitoring Service (CAMS) [23], was implemented in HARMONIE–AROME Cycle 46 [24]. It is not yet used by default. NRT aerosol effects are taken into account both in the radiation and cloud microphysics (Section 3.2.2) parametrizations, while the aerosols represented by the Tegen climatology only influence the radiation parametrization. Details about the coupling of external aerosol data to the HARMONIE-AROME model are given in Appendix A.

The radiation scheme uses aerosol optical depths (AODs) of land (organic matter and sulfates), sea, desert, and urban tropospheric aerosols at a wavelength of 550 nm (AOD550), either from the Tegen climatology or derived from the CAMS NRT data. To obtain three-dimensional AOD550 fields from the NRT aerosol mass mixing ratio (MMR) fields, the mass extinction coefficients suggested by [25] are used. In particular, the mass extinction coefficients representing a relative humidity of 80% are applied. The values of the single-scattering albedo, asymmetry factor, and the spectral distributions of all aerosol inherent optical properties are as previously prescribed within the IFS radiation scheme.

The impact of aerosols on SW radiation was found to be large in the case of desert dust intrusions, as shown in Figure 2 and Figure 3. Generally, in clear sky cases, the SW radiation calculated when NRT aerosols are used was found to be slightly lower compared to that calculated with the Tegen climatology. Figure 2 shows a mostly clear-sky dust example that occurred on the 20 of February 2023 over the Iberian Peninsula. On that day, Saharan dust reduced the global SW radiation by nearly 100 Wm⁻², as observed at some stations. The daily cycle of the areal mean global radiation in the experiments using Tegen and CAMS NRT aerosols is shown in Figure 3. Here, hourly global radiation is compared with observations. The maximum global radiation is overestimated by the experiment with the Tegen aerosol climatology, while in the experiment with NRT aerosols, it is lower and closer to that observed.

Radiation schemes are known to be sensitive to the definition of liquid and ice particle sizes. For a cloud with a given total water content, the optical density is significantly larger when cloud water is distributed among a multitude of small droplets, compared with fewer larger droplets. For the calculation of the cloud liquid droplet effective radius in the radiation scheme of HARMONIE–AROME Cycle 40, constant values of cloud droplet number concentrations (CDNC) were used, with different values over sea and land. CDNC values are also used in the microphysics scheme and described in more detail in Section 3.2.1. In HARMONIE–AROME Cycle 40, different values of CDNC were used in the radiation and microphysics schemes, which was an inconsistency. In HARMONIE–AROME Cycle 46, the radiation and microphysics schemes use the same prescribed base CDNC value with assumed vertical distributions. Optionally, the use of the NRT aerosol data allows us to obtain the 3D CDNC fields and to apply these consistently to both the radiation and cloud microphysics parametrizations (see [24] and Section 3.2.2).

The contributions to HARMONIE–AROME Cycle 46 described in the rest of this section concern LW radiative properties of liquid clouds. HARMONIE-AROME suffers from long-standing issues with forecasting fog: fog is too dense, too cold, and too widespread. As an attempt to advance the solving of this problem, Tosca Kettler performed a study for her Masters thesis in collaboration with KNMI [26]. She investigated the strong cooling of the fog layer and a possible link with LW cloud properties, primarily using MUSC. She suggested improving the calculation of the LW effective emissivity. An equation of the following form is used in HARMONIE-AROME Cycle 46 for the LW effective emissivity $\epsilon$ [27]:

$$\epsilon =1-e^{-0.144(1.2-0.006r_e){\displaystyle \frac{L}{1.66}}},$$

(1)

where L is the integrated cloud liquid water path (LWP) in kg/m², and $r_e$ is the cloud droplet effective radius. The coefficients are used to describe the following: 1.66 is the diffusivity factor [28], 1.2 and −0.006 are used to describe the linear part of the dependency on $r_e$, and −0.144 is an empirically fitted coefficient. In [26], observations from the Cabauw tower and surface instruments were used to re-tune the −0.144 coefficient (the other coefficients were not studied). Figure 4 shows the dependency of the effective LW emissivity on LWP derived from observations and the fit to Equation (1) for other values in addition to −0.144. The effective LW emissivity was derived from the observed downward LW radiative fluxes at the top of the Cabauw tower and at the surface. The LWP was calculated from visibility measurements at six levels in the tower. The −0.144 coefficient was found to be too large Kettler [26] suggested replacing it with the best fit value of −0.096.

As a further improvement of the cloud LW radiation properties calculation, K.P. Nielsen developed a new, spectral parametrization (K. P. Nielsen, personal communication). In this approach, the theoretical Mie-Debye computations for 16 LW spectral bands [14] were used. The core of the new parametrization is given in Equation (2):

$$n_{\lambda }\left(r_e\right)=a_{\lambda }+b_{\lambda }{r}_e+c_{\lambda }{r}_e^2+d_{\lambda }{r}_e^3+{\displaystyle \frac{e_{\lambda }}{r_e}}+{\displaystyle \frac{f_{\lambda }}{r_e^2}}+{\displaystyle \frac{g_{\lambda }}{r_e^3}},r_e>4 \mathrm{\mu }\mathrm{m},$$

(2)

where $n_{\lambda }\left(r_e\right)$ is the spectral band cloud liquid droplet mass absorption coefficient related to the wavelength $\lambda$, and $a_{\lambda }$ - $g_{\lambda }$ are the coefficients for each LW spectral band, as given in

Table 1. LW spectral band coefficients for the Nielsen LW liquid cloud optical property scheme.

LW Band (cm−1)	${\mathbf{a}}_{\mathit{\lambda }}$	${\mathrm{b}}_{\mathit{\lambda }}$	${\mathrm{c}}_{\mathit{\lambda }}$	${\mathrm{d}}_{\mathit{\lambda }}$	${\mathrm{e}}_{\mathit{\lambda }}$	${\mathrm{f}}_{\mathit{\lambda }}$	${\mathrm{g}}_{\mathit{\lambda }}$
1 (10–350)	$-8.521\times {10}^{-3}$	$9.341\times {10}^{-4}$	$-2.611\times {10}^{-5}$	$2.470\times {10}^{-7}$	$4.940\times {10}^{-1}$	−1.496	1.495
2 (350–500)	$-8.737\times {10}^{-4}$	$2.558\times {10}^{-4}$	$-2.210\times {10}^{-6}$	$-2.190\times {10}^{-8}$	$3.678\times {10}^{-1}$	−1.879	3.198
3 (500–630)	$5.902\times {10}^{-2}$	$-2.827\times {10}^{-3}$	$7.173\times {10}^{-5}$	$-6.961\times {10}^{-7}$	$-1.016\times {10}^{-2}$	$-2.569\times {10}^{-1}$	$5.617\times {10}^{-1}$
4 (630–700)	$3.365\times {10}^{-2}$	$-1.835\times {10}^{-3}$	$4.930\times {10}^{-5}$	$-4.967\times {10}^{-7}$	$4.807\times {10}^{-1}$	−2.564	4.461
5 (700–820)	$9.855\times {10}^{-2}$	$-4.580\times {10}^{-3}$	$1.062\times {10}^{-4}$	$-9.421\times {10}^{-7}$	$-3.690\times {10}^{-1}$	1.214	−1.617
6 (820–980)	$3.752\times {10}^{-2}$	$-2.601\times {10}^{-3}$	$6.810\times {10}^{-5}$	$-6.291\times {10}^{-7}$	$7.428\times {10}^{-1}$	−2.333	2.869
7 (980–1080)	$1.204\times {10}^{-1}$	$-5.852\times {10}^{-3}$	$1.278\times {10}^{-4}$	$-1.045\times {10}^{-6}$	$-1.854\times {10}^{-1}$	$5.470\times {10}^{-1}$	$-8.288\times {10}^{-1}$
8 (1080–1180)	$9.444\times {10}^{-2}$	$-3.925\times {10}^{-3}$	$7.286\times {10}^{-5}$	$-5.088\times {10}^{-7}$	$-7.120\times {10}^{-2}$	$1.225\times {10}^{-1}$	$-1.504\times {10}^{-1}$
9 (1180–1390)	$7.449\times {10}^{-2}$	$-2.699\times {10}^{-3}$	$4.307\times {10}^{-5}$	$-2.586\times {10}^{-7}$	$6.176\times {10}^{-2}$	$-3.441\times {10}^{-1}$	$4.269\times {10}^{-1}$
10 (1390–1480)	$7.749\times {10}^{-2}$	$-2.928\times {10}^{-3}$	$5.146\times {10}^{-5}$	$-3.591\times {10}^{-7}$	$4.441\times {10}^{-2}$	$-1.767\times {10}^{-1}$	$6.093\times {10}^{-2}$
11 (1480–1800)	$2.004\times {10}^{-2}$	$-1.134\times {10}^{-3}$	$2.468\times {10}^{-5}$	$-1.993\times {10}^{-7}$	$8.347\times {10}^{-1}$	−2.580	3.197
12 (1800–r2080)	$-5.067\times {10}^{-2}$	$1.808\times {10}^{-3}$	$-3.382\times {10}^{-5}$	$2.552\times {10}^{-7}$	1.515	−2.572	1.520
13 (2080–2250)	$-5.707\times {10}^{-2}$	$2.044\times {10}^{-3}$	$-3.788\times {10}^{-5}$	$2.816\times {10}^{-7}$	1.575	−1.810	$-4.621\times {10}^{-1}$
14 (2250–2380)	$-7.484\times {10}^{-2}$	$2.538\times {10}^{-3}$	$-4.423\times {10}^{-5}$	$3.096\times {10}^{-7}$	1.882	−3.059	$4.880\times {10}^{-1}$
15 (2380–2600)	$-7.583\times {10}^{-2}$	$1.924\times {10}^{-3}$	$-2.835\times {10}^{-5}$	$1.873\times {10}^{-7}$	2.461	$-1.015\times {10}^1$	$1.415\times {10}^1$
16 (2600–3250)	$7.079\times {10}^{-2}$	$-6.443\times {10}^{-4}$	$-1.294\times {10}^{-5}$	$1.528\times {10}^{-7}$	$-2.350\times {10}^{-1}$	$6.213\times {10}^{-1}$	$-5.661\times {10}^{-1}$

. The calculation of the LW cloud emissivity from the spectral band mass absorption coefficients by Equation (2) confirmed that it was indeed too high in HARMONIE-AROME.

Spectral absorptance calculations are shown in Figure 5. The corresponding values calculated from the emissivity by the Smith and Shi [27] parametrization (purple continuous line) and the Kettler scheme (green continuous line [26]) are presented. It can be seen that the Kettler parametrization [26] fits the results of the new theoretically based method well for wavelength bands 8, 9, and 10, which represent the smallest, most common cloud droplet sizes. These spectral bands are the most important for the cloud LW radiative effects, because they are least affected by greenhouse gas absorption. However, in an atmosphere with highly variable water vapor concentrations, other wavelength bands are also important. The new scheme improves the representation of the LW cloud properties compared with the scheme based on the non-spectral LW cloud emissivity parametrization. It is a default now in HARMONIE-AROME Cycle 46.

3.2. Cloud Microphysics

The core of the parametrization of cloud microphysics in HARMONIE-AROME is a one-moment bulk scheme containing three classes of ice, based on developments originally carried out in Meso-NH [29,30]. This ice classification is commonly referred to as “ICE3”. The three classes of ice (solid hydrometeors) are snow, cloud ice, and a combination of hail and graupel. For liquid and gaseous hydrometeors, three classes are taken into account: rain, cloud water, and water vapor. All of these hydrometeor classes are represented by prognostic mixing ratios and advected by the model dynamics; horizontally by the semi-Lagrangian scheme [3] and vertically by the sedimentation process [31]. The grid box-wise particle sizes are estimated from a generalized gamma distribution.

This section contains four subsections covering the developments of (1) a better representation of cloud droplet number concentration that affects fog and clouds; (2) using NRT aerosol data for cloud microphysics calculations; (3) the OCND2 option to improve forecasts in winter, especially for cases with a stable boundary layer; and (4) the ICE-T scheme, to improve the representation of supercooled liquid water.

3.2.1. Cloud Droplet Number Concentration

In the current default version of HARMONIE-AROME Cycle 46 (where NRT aerosols are not used), microphysically important variables are prescribed. The values are either height dependent, as described in Contreras Osorio et al. [32], or approximated process-wise. One such important meteorological variable is the cloud droplet number concentration (CDNC). It is used in the parametrization of various processes leading to the growth of activated cloud droplets to liquid and solid precipitation, see further details in Section 3.2.3 and Section 3.2.4.

In HARMONIE-AROME Cycle 43, constant values of CDNC were used, with a value of 500 cm⁻³ for urban areas, 300 cm⁻³ over land, and 100 cm⁻³ over the sea/ocean. With these values, fog and cloud condensate in the lowest thickest clouds were over-estimated. An example of the problems with forecasting low thick clouds is shown in Figure 6. In this figure, a meteosat second generation (MSG) visible satellite image is shown, along with the MSG Seviri cloud water path product from the Dutch Meteorological Service, KNMI, and cloud water condensate fields from two HARMONIE-AROME Cycle 43 experiments, for 12 UTC on 8 July 2019. One experiment used the default configuration in HARMONIE-AROME Cycle 43, and the other one used CDNC values of 50 cm⁻³ everywhere, as well as the LW effective emissivity suggested by Kettler [26]. One can see from Figure 6 that the integrated cloud water was overestimated in the default HARMONIE-AROME Cycle 43. A CDNC value of 50 cm⁻³, with a new LW effective emissivity, gave much better results.

In HARMONIE-AROME Cycle 46 CDNC has a vertical dependence with height, with the same profile over land and sea. The reason for this is to eliminate the artificial reduction in stratiform precipitation that occurred in HARMONIE-AROME Cycle 43 at the land sea boundaries, due to the use of different CDNC values. At the lowest model level, CDNC is multiplied by a factor of 0.25. This is to improve predictions of low visibility and fog. The current configuration of HARMONIE-AROME Cycle 46 shows improvements in the cloud water prediction.

3.2.2. Effects of CDNC Obtained from the NRT Aerosol Fields

The use of NRT aerosol data in HARMONIE-AROME (as described in Appendix A, see also [24] for more details) allows us to calculate the CDNC instead of prescribing it. The estimation of CDNC from the aerosol MMR data is based on the Köhler theory. Hydrophilic aerosols (sea salt, sulfates, nitrates, ammonium, and hydrophilic organic matter and black carbon) are activated under supersaturated conditions. The supersaturation within clouds is calculated based on thermodynamic variables and the vertical velocity. Case studies [24,33] show that the precipitation amount, and in particular, its phase (snow, graupel, or rain), is sensitive to the use of aerosol MMRs in the cloud microphysics.

Since cloud distribution affects surface radiation fluxes, we can use observations of SW radiation for an indirect validation of cloud parametrizations. Here, we use observations of SW radiation from 20 synoptic sites in Ireland to validate the impact of the use of NRT aerosols on cloudiness. Figure 7 shows distributions of the observed and modeled clear sky index (CSI) for experiments with NRT aerosols (CAMSNRT) and without NRT aerosols (i.e., with prescribed CDNC, Tegen experiments). CSI is the ratio of the observed or modeled global SW radiation to the clear sky global SW radiation (the lowest values thus correspond to cloudy conditions, while values close to 1 represent clear skies). The results for a 2 week summer period (1 June to 14 June 2018) are shown in Figure 7a and for a winter period (3 February to 17 February 2020) in (b). A clear overestimation of the low CSI values in the experiment with prescribed CDNC can be seen in both seasons. This is consistent with the overestimation of thick clouds mentioned in Section 3.2.1 and shown in Figure 6. When NRT aerosols are used, this overestimation of low CSI values is not present. Figure 8 shows the global SW radiation bias distributions for the two experiments. The positive and negative biases are plotted on the positive axis here, to highlight whether the experiments result in more positive or negative biases overall. It can be noticed that in the experiment with prescribed CDNC, there are negative biases overall, while when NRT aerosols are used, the biases are overall positive. Hence, the NRT aerosol usage results in a general overestimation of the SW radiation. It is important to note here that in these experiments, the impact of the NRT aerosols on both the radiation and the microphysics schemes is seen.

3.2.3. OCND2: Clouds in Cold Conditions

As described in [1,34], the ICE3 scheme was supplemented with an option called OCND2 to improve simulations in winter, especially during stable boundary layer conditions over Northern Europe. Under these conditions, clouds contain supercooled liquid rather than ice. This situation is difficult for the model to reproduce. Moreover, sometimes the model suffers from spurious fog at temperatures below −20 °C. Here, we list the aspects of the ICE3 microphysics parametrization modified under the OCND2 option in order to address weaknesses in HARMONIE-AROME Cycle 46:

• The separation between liquid water and ice water processes was improved. This means that the statistical cloud scheme (See Section 3.3.3) only deals with cloud liquid water, including cases when temperatures are below freezing. All ice processes are described by the OCND2 version of the ICE3 scheme.
• Evaporation/deposition of cloud ice is a conversion between ice and vapor and not between ice and liquid.
• The deposition rate of the ice species was reduced.
• The cloud cover, from the point of view of the forecast users (the public), was modified to account for the lower optical thickness of ice clouds compared to water clouds.
• The ice nucleating particle number concentration is reduced when the temperature is between 0 °C and 25 °C. The main purpose of this is to slow down the conversion from cloud liquid water to ice, snow, or graupel.
• To support the production of supercooled rain, threshold values were introduced for converting supercooled rain into graupel or snow.
• Calculations of saturation pressure are avoided when the saturation pressure is near or above atmospheric pressure. This is implemented for computational reasons and affects calculations in the stratosphere only.
• A bugfix was implemented to solve a technical problem, which revealed spurious ice-clouds (see Figure 9). The ICE3 scheme should be active when the amount of any non-vapor water species exceeds some low threshold or when the air temperature is below freezing. Unfortunately, for technical reasons, this did not always happen when the second criteria was satisfied. Due to that, areas with sufficient water species were surrounded by areas with little cloud ice water, as shown in Figure 9.

With the OCND2 changes, the microphysics scheme is more “liquid friendly”, which means that supercooled droplets freeze more slowly. However, it is still not “liquid friendly” enough. There are cases where supercooled clouds seem to disappear too quickly. A further implementation, aimed at improving the representation of supercooled liquid water in HARMONIE-AROME, is the ICE-T scheme, described in the next subsection.

3.2.4. ICE-T: To Improve the Representation of Supercooled Liquid Water

The development of ICE-T was initiated during a study of ice loads on transmission lines from atmospheric icing. The study’s hypothesis was that HARMONIE-AROME (even with OCND2) glaciates the clouds prematurely, since the calculated loads on transmission lines, when applying the supercooled liquid water (SCL) from the model, was underestimated [35]. Elements from the Thompson microphysics scheme [36] were implemented in the ICE3 microphysics in addition to the OCND2 changes. The Thompson microphysics scheme [36] is also used in the Weather Research and Forecasting (WRF) model. The name of the new scheme is ICE-T, which is the combination of ICE3 and the Thompson scheme. The scheme is described in detail in [37]. Its most prominent features are listed below:

• Stricter conditions for ice nucleation.
• Less efficient collision–collection of liquid water by solid hydrometeors.
• A variable rain intercept parameter, which allows for smaller droplets when condensation and coalescence are the primary sources.

ICE-T was tested in two studies regarding atmospheric icing on power lines [35] and aircraft [38]. The results show a clear shift towards more SCL water. Both studies found that the increased SCL water was in better correspondence with observations of ice loads, measurements of atmospheric content of liquid and ice water by satellites, and pilot reports of aircraft icing. However, the atmospheric content of ice species also increases. The reason is a shift from graupel to snow. Since snow has a slower terminal fall speed compared with graupel, the residence time is increased, and hence, the accumulated amount is increased. There is an increase in surface precipitation in the form of snow and a decrease in the form of graupel. The total precipitation is reduced by a few percent, and the precipitation pattern is shifted from the upstream to the lee side of topography. The shift in the precipitation pattern will be explored further.

Recently, ICE-T was validated using data from a helicopter icing measurement campaign launched from Alta, Norway in April 2023. The campaign was led by Airbus in order to test their helicopters’ ability to fly through heavy icing conditions. On 19 April 2023, heavy icing conditions inside altocumulus lenticularis clouds occurred over mountainous areas in the vicinity of Alta. The helicopter was carrying a Cloud Droplet Probe for measuring hydrometeors in the range 2–50 µm, which are essentially cloud droplets. During an afternoon flight they measured liquid water content of mostly 0.8 gm⁻³ and above. The highest value measured was about 1.3 gm⁻³. Two parallel HARMONIE-AROME Cycle 46 simulations were performed for the validation of the ICE-T scheme against these observations: one with the default ICE3 scheme with the OCND2 option (reference experiment) and the other with the ICE-T scheme. The simulations were initialized on 19 April 00 UTC, with no upper-air data assimilation. Figure 10 shows the cloud water content (gm⁻³) for both simulations and the difference between these at the level of approximately 820 hPa, where the helicopter encountered the heaviest icing conditions. The figure shows that the ICE-T simulation gives higher cloud water content compared to ICE3, and this is also closer to what was observed from the helicopter.

3.3. Shallow Convection, Turbulence, and Statistical Cloud Scheme

A comprehensive integral revision of three parametrization schemes in HARMONIE-AROME was carried out by [39]: the statistical cloud scheme, the turbulence scheme, and the shallow cumulus convection scheme. These schemes collectively parametrize boundary layer processes. Hence, in HARMONIE-AROME Cycle 46, the cloud, turbulence, and convection schemes differ significantly from the implementation described in [1]. A response to these updates on subgrid- and grid-scales for a cold air outbreak case is described in [40]. We describe these updates briefly in the subsections below. A further implementation in the turbulence scheme, which concerns wind farms and is not related to this revision, is described separately in Section 3.3.4.

3.3.1. Shallow Convection

At a grid spacing of 2.5 km, deep convection is expected to be roughly resolved and explicitly represented by the model’s non-hydrostatic dynamics. HARMONIE–AROME, therefore, does not parametrize deep convection. However, shallow convection still needs to be parametrized. The shallow convection scheme applies a dual mass flux framework [41], in which two updrafts are distinguished: a dry updraft that does not enter the cloud layer and a moist updraft that reaches the lifting condensation level and continues its ascent in the cloud layer. The updated shallow convection scheme is described in detail by [39], where all the modifications since the HARMONIE-AROME Cycle 40 configuration [1] are included. One of the most important modifications concerns the coupling of the convection and the turbulence schemes via the so-called energy cascade term. The lateral mixing term, from the prognostic mass flux vertical velocity variance equation [42], is used as a source term in the turbulent kinetic energy (TKE) equation. This particularly enhances the subcloud-to-cloud layer transport. This is in better correspondence with large eddy simulation (LES) results for shallow convection, as shown in [39] and illustrated by Figure 11. This figure shows the total turbulent transport of moisture ($\overline{w^{\prime }{r}_{t,tot}^{\prime }}$) according to the Dutch Atmospheric Large-Eddy Simulation model (DALES) [43], HARMONIE-AROME with an energy cascade, and HARMONIE-AROME without an energy cascade. In the case of HARMONIE-AROME, $\overline{w^{\prime }{r}_{t,tot}^{\prime }}$ is the sum of the transport by the turbulence and convection schemes, whereas for the DALES model, the total turbulent transport consists of the resolved turbulent transport and a small sub-grid part. These total turbulent fluxes are most important because the vertical divergence of them determines the tendencies of the corresponding prognostic variables. Figure 11 shows that like the DALES model, the version of HARMONIE-AROME with the energy cascade term describes the increasing moisture transport with height (i.e., the drying) in the sub-cloud layer. As a result, the ventilation of the sub-cloud layer and moistening of the cloud layer are represented much better. Simulations were performed for the Atmospheric Radiation Measurement (ARM) shallow cumulus case [44].

3.3.2. Turbulence

In this section, we describe changes to the turbulence scheme especially aimed at improving the forecasting of low clouds. Starting with HARMONIE-AROME Cycle 36, the HArmonie with RAcmo TUrbulence (HARATU) [1,39,45] scheme replaced the Cuxart–Bougeault–Redelsperger (CBR) scheme [46] for describing turbulence. Both schemes apply a prognostic equation for TKE and a diagnostic length scale. Although HARATU significantly improved several aspects of turbulence (see [1,47,48]), it also contributed to an underestimation of low cloud cover and an overestimation of cloud base heights. In the comprehensive integral approach described in [39], the turbulence scheme was revised to improve the forecasts of low clouds in particular (see Figure 12). This figure shows the frequency bias of cloud base height classes for one month of simulations using HARMONIE-AROME Cycle 40 [1] (a) and with the modifications of [39] (b). The blue, green, and orange lines refer to +3, +24, and +48 h forecasts, respectively. In the reference HARMONIE-AROME Cycle 40, for +24 h forecasts of the cloud base height around 178 feet (x-axis), the fraction of predicted cases relative to the number of the observed cases is less than 20% (y-axis). For the updated model version, shown in (b), the prediction of events with a low cloud base height has clearly improved. Two modifications to the turbulence scheme are most relevant to the improved prediction of low clouds. The first one is related to the re-tuning of the scheme based on the Monin–Obukhov similarity theory, following [49,50]. The second one concerns changing the free asymptotic mixing length. One of the key consequences of these modifications is a better conservation of atmospheric inversion strengths. Consequently, stratocumulus clouds are better preserved (and the triggering of deep precipitating convection is also improved) [39].

3.3.3. Statistical Cloud Scheme

In HARMONIE-AROME, liquid water and mixed-phase clouds, which commonly appear as low and medium clouds, are described by a statistical cloud scheme. Ice clouds, which are mainly high clouds, are parametrized in a relative humidity scheme, see Section 3.2.3. In the statistical cloud scheme, the cloud fraction and liquid water content can be derived from the sub-grid variability of temperature and moisture [51]. A detailed description of the statistical cloud scheme, including several fixes and modifications to the Cycle 40 version of HARMONIE-AROME, can be found in [39]. The main modifications concern the following:

• The derivation of the thermodynamic coefficients. This is performed in a proper way now.
• Adding of the covariance term for temperature and humidity.
• The dissipation length scale in the cloud scheme is now consistent with the one in the turbulence scheme.
• The description of the dissipation of the variances. Now it is more consistent with the literature.
• Removing of the erroneous factor 2 for the convective contribution to the variance.

Although the updated statistical cloud scheme is based on a more sound physical basis and several errors are removed, results for the idealized ARM cumulus case [44] resulted in an even more pronounced underestimation of the variance, see Figure 13. The variance is especially underestimated in the upper part of the cloud layer, leading to an even greater underestimation of cloud fractions at these heights [39]. Note also that for the time stamp shown in the figure, the convective fluxes, which are an important source of s^′2 in the LES and in HARMONIE-AROME, match closely (see [39]). To solve this problem, we are currently investigating a physically plausible approach for introducing the dependency of the dissipation time scale of convection on height.

3.3.4. Wind Farm Parametrization

As investment in renewable energy is increasing, wind turbine farms cover large areas and increasingly have become a part of the European landscape. Large wind turbines affect atmospheric motions, inducing additional turbulence, typically in layers higher than the surface layer but still in the boundary layer. In order to model this effect, a wind farm parametrization (WFP) was implemented in HARMONIE-AROME Cycle 46. This parametrization alters the tendencies of momentum U and TKE by

$${\displaystyle \frac{\partial U_k}{\partial t}}=-{\displaystyle \frac{1}{2}}{C}_{\mathrm{T}}{U}_k\left|\overrightarrow{V_k}\right|A_k{\Delta }_k^{-1},$$

(3)

$${\displaystyle \frac{\partial {\mathrm{TKE}}_{\mathrm{k}}}{\partial t}}={\displaystyle \frac{1}{2}}{C}_{\mathrm{TKE}}{\left|\overrightarrow{V_k}\right|}^3{A}_k{\Delta }_k^{-1}.$$

(4)

These equations are solved for each model level k. Here, U is the two wind components (u and v), $|\overrightarrow{V}|=\sqrt{u^2+v^2}$; A is the rotor area; ${\Delta }_k$ is the volume of the grid cell, $C_{\mathrm{TKE}}$ = $C_{\mathrm{T}}-$$C_{\mathrm{P}}$; $C_T$ is the trust coefficient; and $C_P$ is the power coefficient. Trust and power coefficients are defined for each wind speed and wind turbine. The way in which $C_{\mathrm{TKE}}$ is calculated means that the remainder of the energy extracted from momentum, which is not converted to power, is transferred to TKE.

The implementation of this parametrization is described in Van Stratum et al. [52]. There, an evaluation of the impact is shown, using one year of HARMONIE-AROME simulations. An example of the modeled wind speed at hub height in the German Bight is shown in Figure 14. The parametrization produces wakes from the wind turbines with a velocity deficit of up to 2 ms⁻¹, also visible in aircraft observations. The WFP needs additional external data about the location and parameters of the wind turbines, which are usually available locally. Because of that, this parametrization is not a default in HARMONIE-AROME Cycle 46.

4. Surface Physics

In HARMONIE-AROME, the SURFEX modeling platform [54] is used to model the surface of the Earth. This platform is mainly developed and maintained by Météo-France. In HARMONIE-AROME Cycle 46, SURFEXv8.1 (hereafter SURFEX) is used with several technical and scientific modifications applied for use in HARMONIE-AROME. SURFEX contains physical models of varying levels of complexity for four main types of surface (called tiles): natural land surfaces (soil and vegetation), urban areas, inland waters, and oceans. SURFEX includes the surface layer, which, in NWP models, is traditionally considered to be the layer between the lowermost model level and the surface of the Earth. This layer is an interface between the atmospheric and surface parts of the model and provides the exchange of turbulent fluxes. The nature tile can be further divided into sub-types (called patches) to describe different kinds of vegetation, e.g., crops, forestry, and rocks (with no vegetation). SURFEX also includes the software needed to initialize models and to prepare land-cover information, including all of the necessary external parameters such as leaf area index (LAI), albedo, tree heights, and lake depths. For the nature tile, the vegetation parameters are specific to each patch. Fluxes and diagnostic variables are calculated for all tiles and patches and then aggregated for each grid box using their corresponding fractions.

The following subsections outline developments carried out within the HARMONIE-AROME framework in the areas of physiography, urban and nature tiles, snow melt, the stable boundary layer, orographic enhancement of momentum fluxes, and the sea and inland water tiles.

4.1. Physiography

The land cover characteristics in HARMONIE-AROME Cycle 46 are described using the ECOCLIMAP Second Generation (ECOSG) land cover map (https://opensource.umr-cnrm.fr/projects/ecoclimap-sg, accessed on 27 October 2024) developed by Météo-France, which has a primary resolution of ca. 300 m. This is an advancement compared with HARMONIE-AROME Cycle 40, where the 1-km resolution ECOCLIMAPv.2.2 [55] was used. ECOSG was corrected for use in HARMONIE-AROME Cycle 46 by applying the Randolph Glacier Inventory [56] to give a better representation of the glaciers in Norway, Iceland, and Svalbard (they are larger in the original ECOSG). As well as this, land cover updates were suggested over Iceland, Greenland, Svalbard, and the Faroe Islands, using freely available land cover data from different sources [57,58,59]. An example of the improvements is shown in Figure 15.

To represent topography, the GMTED2010 dataset (https://topotools.cr.usgs.gov/gmted_viewer/viewer.htm, accessed on 27 October 2024) is used in HARMONIE-AROME Cycle 46. It has been corrected using data from the ArcticDEM [60] and Copernicus GLO-90 (https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.1, accessed on 27 October 2024) projects. For the soil properties, we rely on the SoilGrids database [61], which has a primary resolution of ca. 300 m. For lake depths, we use the Global Lake Database v.3 [62], with a primary resolution of ca. 1 km, and for the initialization of lake variables during the first model run (called a cold start), we use the global lake climatology by [63].

4.2. The Urban and Nature Tiles

The urban tile is simulated by the Town Energy Balance (TEB) model [64]. Over the nature tile, the Interaction between Soil, Biosphere, and Atmosphere (ISBA) scheme simulates the energy and water exchange between the land surface and the atmosphere and the evolution of temperature and moisture within the soil. Currently, HARMONIE-AROME Cycle 46 uses the force-restore (ISBA-FR; Boone et al. [65]) version of the scheme. A one-layer snow scheme [66] is used to describe the evolution of snow cover. Diagnostic temperature and humidity at 2 m and wind speed at 10 m, over land (as well as over other surface types) are calculated by interpolation in the surface layer, with a dependency on turbulent fluxes. In HARMONIE-AROME Cycle 46 compared with Cycle 40, we increased the number of nature patches from 1 to 2 so that we now have distinct patches for low and high vegetation types.

While using two patches in combination with the high-resolution ECOSG land cover map, we experienced an increase in 10 m wind speed forecasts over areas covered mainly by low vegetation. This resulted in a positive bias in 10 m wind, e.g., over Ireland. Only pure land cover types exist in ECOSG, whereas previous versions of ECOCLIMAP had mixed land cover types. At fine resolution, often only low vegetation is present in a nature tile, which leads to quite a smooth surface. As a result, most natural areas in Ireland, for example, are represented as grassland. In reality, even in low-vegetation areas, there are usually some trees, which affect the airflow. Fields are often surrounded by hedgerows and sporadic trees. The absence of such features in the underlying land cover dataset means there is an absence of roughness elements to slow down near-surface winds in the model forecasts. This effect is not important with the one-patch “mean” vegetation approach, but with ECOSG, and the division of vegetation into low and high categories, this becomes important. The option to parametrize this effect by increasing the vegetation height is included in HARMONIE-AROME Cycle 46 (a so-called LFAKETREE option): 10% of the low vegetation (grassland and crops) in each grid box can be replaced by trees of height 10 m. This option only affects the average roughness length of the low vegetation patch through the logarithmic averaging of roughness length and does not affect the rest of ISBA. Figure 16 shows verification scores for two test periods, with different physics options: HARMONIE-AROME Cycle 40 (red), HARMONIE-AROME Cycle 43 default (green), and HARMONIE-AROME Cycle 43 with “LFAKETREE” (blue). Only results from the 00 UTC forecast cycles are shown in order to highlight any diurnal patterns. From Figure 16, it is evident that the LFAKETREE option reduces the wind-speed bias.

For forestry, the surface roughness is calculated as a function of tree height. ECOSG includes a supplementary tree height dataset. However, at northern latitudes, the tree height values led to roughness lengths that were too high and, hence, to an underestimation of surface wind speeds. To circumvent this problem, a tree height scaling with latitude was introduced.

4.3. Snow Melt Adjustment

The one-layer snow scheme [66] has a problem with melting very small snow amounts. This issue is caused by the use of the snow fraction $f_{sn}$ in the snow melt calculation, which is given as follows:

$$MELT=f_{sn}{\displaystyle \frac{\Delta T}{[C_p\times c_m\times \Delta t]}}$$

(5)

where $MELT$ is the snow melt in kg·m⁻² · s⁻¹, $\Delta T$ is the energy available for melting expressed as a temperature excess over 0 °C, $C_p$ is the snow thermal coefficient, $c_m$ is the melting energy of ice, and $\Delta t$ is the time step. When the snow fraction becomes very small, the snow melt also becomes very small. As a result, the last few millimeters of snow take an unrealistically long amount of time to melt. The impact of this on the surface fluxes is small, as $f_{sn}$ is small when the snow amount is small. However, the snow does not disappear from the model output quickly enough, which is deceiving for forecasters. The suggested solution was to limit $f_{sn}$ to some minimum value. This value was found from experiments with MUSC. A minimum value of 0.25 gives reasonable results, as shown in Figure 17.

4.4. Stable Boundary Layer

Many NWP models experience difficulties in correctly representing the boundary and surface layer under stable conditions, which influences the evolution of 2 m temperature. This is known as the stable boundary layer (SBL) problem [67,68]. Too much heat and momentum transport results in a warm bias in near-surface temperatures and reduced low-level wind speeds [69]. Halting this transport triggers a runaway cooling regime that causes unrealistically cold surface temperatures [70,71]. Turbulent flux calculations in the SBL are frequently the target of model tuning to enhance performance. Such tuning efforts have also continued for HARMONIE-AROME. The main tuning parameter was a threshold for the gradient Richardson number. In the model, there is an upper limit for this threshold called XRIMAX. This leads to a limitation in the calculations of the drag coefficients and, hence, in the turbulent fluxes.

In HARMONIE-AROME Cycle 40, the default value of XRIMAX was 0, which prevented the formation of a stable boundary layer. Experiments with XRIMAX > 0 (e.g., 0.2 or 0.5) showed potential for improved performance in cold stable situations at inland stations in Scandinavia. However, this also resulted in unrealistic temperature drops, particularly at coastal stations in the north and at inland stations in mountainous regions, where the roughness length was low and when the ground was covered with snow. Homleid [72] suggested modifying the dependency of the drag coefficients on the gradient Richardson number. Her suggestion involved restricting the drag coefficients to be equal to neutral values, for a certain interval of stability characterized by the gradient Richardson number. This resulted in an improvement in the scores in some areas. In general, model performance is very sensitive to the different options used for calculating the drag coefficients and to the values of XRIMAX. However, it is impossible to find an optimal solution that would result in a performance that is equally good for all geographical conditions.

Algorithms used to interpolate meteorological variables to near-surface values can differ, which also affects the verification scores. Near-surface scores should not be the only criteria used for making decisions; the model performance at levels higher than the surface layer is also very important. A detailed depiction of the impacts of different options on stable boundary layer behavior can be found in Kähnert et al. [73]. The default value of XRIMAX in the reference version of HARMONIE-AROME Cycle 46 is 0.2, but this value does not give the best scores in domains dominated by snow in winter.

4.5. Orographic Enhancement of Momentum Fluxes

The effect of small-scale orographic variability on turbulent momentum fluxes over a complex terrain can be taken into account in SURFEX by using a choice of options. One option uses the orographic roughness concept [74], possibly with directional components based on Georgelin et al. [75]. The problem with the orographic roughness approach is that the roughness length can reach large values, even higher than the thickness of the surface layer, which is assumed to extend from the surface to the lowest model level (currently ca. 14 m). This may lead to nonphysical behavior of the surface drag (a decrease) and low-level winds (an increase) over rough terrain. Smaller, but still too large, roughness values tend to lead to a drag that is too large, which retards the near-surface winds excessively. An alternative option is a parametrization of the non-separated sheltering effect due to airflow over hills and mountains [76] in the form suggested by Beljaars et al. [77].

A third option, called OROTUR, was suggested and implemented according to [78]. This is another realization of the Wood et al. [76] idea. The OROTUR scheme calculates the orographic stress from the turbulent momentum stress proportional to the subgrid-scale orography variance. The proportionality coefficient is, in turn, calculated from the wind speed at the lowest model level, scaled with some constant wind speed and the grid box size $\Delta x^2$. The idea behind the wind scaling is to increase the drag for the weakest winds by accounting for the surface layer wind shear. Scaling with $\Delta x^2$ allows the orographic variations to be roughly related to the steepness of the subgrid-scale slopes in each grid box. The grid-scale orographic stress is then added to the grid-scale turbulent stress so that no tiles or patches are considered. The method includes two tunable constants. The default version of HARMONIE-AROME Cycle 46 does not use any of the small-scale orography parametrizations. However, in the Spanish Meteorological Agency (AEMET), OROTUR is applied operationally. Calvo et al. [79] reported that the scheme contributed to a reduction in the positive bias in 10 m mean wind speeds and gusts. In HARMONIE-AROME, vertically propagating buoyancy waves from the orography are not parametrized because their generation is believed to be described by the fine-resolution non-hydrostatic dynamics of the model. Dissipation of the waves is assumed to be accounted for by the turbulence parametrization [80].

4.6. The Sea Tile

To represent the sea tile, the sea surface temperature (SST) is initialized from external sources and does not change during the forecast. External fields usually come from ECMWF’s IFS, which, in turn, uses the OSTIA analysis (https://ghrsst-pp.metoffice.gov.uk/ostia-website/index.html, accessed on 27 October 2024) with some modifications. Other SST sources are also used; for example, MetCoOp uses results from the NEMO ocean model [81] running at the Swedish Meteorological and Hydrological Institute (SMHI) over the Baltic Sea (and over the large Swedish lakes of Vänern and Vättern, which are treated as sea water). The sea ice concentration (SIC) is also kept constant and initialized from external sources, usually IFS. Technically, the system allows the SST and SIC fields to be updated from external sources during the integration of the atmospheric model, though this functionality is not currently used in NWP-mode. However, this is applied when running HARMONIE-AROME as a climate model; see, for example, Belušić et al. [82].

The Exchange Coefficients from the Unified Multicampaign Estimates (ECUME) scheme by [83] is used to calculate turbulent fluxes in the surface layer over unfrozen seas. There are two versions of this scheme available in SURFEX, namely, ECUME and ECUME6 [54]. They were developed using observations from different measurement campaigns and use slightly different parameters and formulations. Extensive testing was performed in order to assess their impacts, with observations from the EUREC⁴A project (https://eurec4a.eu/, accessed on 27 October 2024) used in the evaluation process. Comparisons with the EUREC⁴A observations were performed in phase-space; see [84] for the explanations. Figure 18 shows some results from this study, where the ECUME6 scheme overestimates the latent heat flux. However, the bias is small for the ECUME scheme. These results were also confirmed by climate run validations (Oskar Landgren and Bert van Ulft, personal communication). The default scheme in the reference version of HARMONIE-AROME Cycle 46 is, therefore, ECUME.

Similarly to HARMONIE-AROME Cycle 40 [1], sea ice is represented by a one-dimensional thermodynamic parametrization scheme called SICE [85]. However, the default configuration of this scheme has been updated based on the findings of [86] to include prognostic formulations for the ice thickness. As an option, SICE allows the evolution of snow on top of the sea ice to be represented explicitly by means of the ISBA-ES snow scheme. In HARMONIE-AROME Cycle 46, SICE uses an updated version of ISBA-ES available in SURFEX, with 12 snow layers (compared to the 3 layers in HARMONIE-AROME Cycle 40).

4.7. The Inland Water Tile

In HARMONIE-AROME Cycle 40, the surface temperature of inland water (the temperature of lakes, reservoirs and rivers) was initialized using the deep soil temperature and kept constant during the forecast; see [1]. HARMONIE-AROME Cycle 46 uses the FLake [87] lake model to represent the inland water tile. FLake is parametric: it uses a self-similarity concept and predicts the characteristics of the temperature profile in a water column, in ice, and in snow. This approach incorporates much of the essential physics and offers a good compromise between physical realism and computational cost. The temperature profile in a water column is represented by the mixed layer and the thermocline. The mixed layer depth, the mean water temperature, and the bottom temperature of the lake are predicted. The integral of the temperature profile in the thermocline (the so-called shape factor) is also predicted. The surface temperature of the water is diagnostic, and the temperature profiles in snow and in ice are assumed to be linear. The model predicts the depths and temperatures of snow and ice. The snow block uses modifications by Semmler et al. [88]. FLake contains a bottom sediments block, which is switched off by default in HARMONIE-AROME. The model is designed to be applied to fresh and brackish water bodies, although SURFEX uses it to simulate all lakes (and assumes a certain error for saline lakes).

FLake results in an improvement in the representation of the inland water tile, compared with the previous implementation, as illustrated by the example in Figure 19. In this figure, scatter-plots of the simulated and observed values of 2 m temperature and specific humidity are shown, for two HARMONIE-AROME experiments for spring 2016, with observations from three stations in the vicinity of Lake Ladoga. For both 2 m temperature (plots (a) and (b)) and specific humidity (plots (c) and (d)), more points lie along the diagonal when FLake is used.

5. Upcoming Developments in HARMONIE-AROME

While Cycle 46 of HARMONIE-AROME will soon be used operationally, developments will continue in Cycle 49 and beyond. A flavor of these developments is provided in the following subsections in the areas of dynamics and model configuration, radiation, cloud microphysics, shallow convection, and surface physics.

5.1. Dynamics and Model Configuration

The horizontal resolution of NWP models has steadily increased over recent decades. Global models now approach the kilometer scale [89]. Correspondingly, the resolution of limited-area models will continue towards sub-kilometer scales, although this effort faces numerous challenges [90,91].

Model performance at increasing resolution has been an active area of research and development within the ACCORD consortium in recent years. A number of countries now are running experimental sub-kilometer forecasting suites using HARMONIE-AROME. Achieving numerical stability through the choice of an appropriate dynamics configuration is one of the first concerns, with a number of possible avenues available. Stability will be improved by (i) the use of a quadratic or cubic spectral grid, as discussed earlier, due to a reduction in the formal spectral resolution; (ii) off-centering of nonlinear terms in the time discretization; and (iii) in particular, the second-order averaging introduced by [92].

More fundamental considerations of the time discretization scheme may be necessary at extremely high resolutions. As outlined in Section 2, the general ICI approach is designed to improve stability. More iterations will allow for stable integrations but at increasing costs. So far, preliminary testing with the single-step SETTLS-SI approach of HARMONIE-AROME at grid resolutions of around 500 m was successful. While stable forecasts were obtained below this resolution, in the case of complex orography and extreme winds, this sometimes required undesirable parameter choices. For example, extremely cold values are used for the ’reference elastic temperature’ in the SI reference state; see [3,93] for further discussion of this. More investigation is required, and within the wider ACCORD community, considerable effort is underway to address these challenges; for example, see [94,95].

The increase in the horizontal resolution desirably has to be followed by a corresponding increase in the vertical resolution. Therefore, two options for a 90-level vertical grid in HARMONIE-AROME were developed and are being tested (Figure 20c,d). The first option is based on the operational AROME-France 90-level vertical grid (here referred to as $\mathrm{MF}\_90$), with the highest model level at 10 hPa (around 31 km altitude) in the stratosphere and the lowest level around 5 m above the Earth’s surface (with a terrain-following vertical coordinate). MetCoOp developers suggested another option by slightly modifying the $\mathrm{MF}\_90$ vertical grid to place the lowest model level 10 m above the surface (this option for the vertical grid is referred to as $\mathrm{MC}\_90$). The denser vertical grid, especially in the part of the troposphere closest to the surface, improves the representation of the boundary-layer dynamics, thermodynamics, and transport of momentum, energy, and matter by resolving more smaller-scale processes. This, in turn, has a consequence with low-level clouds and fog forecasts being more realistic. An example of the effect of the vertical resolution on cloud fraction is shown in Figure 20a,b. In this figure, vertical cross-sections of cloud fraction (a) along the blue dashed line drawn in a satellite image (b) are shown. The satellite image (b) provides information about the observed cloudiness, showing boundary-layer horizontal convective rolls. The figure shows HARMONIE-AROME simulations of cloud fraction for vertical grids with 65 and 90 (MC_90) levels. The simulation using the grid with 90 vertical levels gives less cloudiness (but more organized structures) than the grid with 65 vertical levels, which seems more realistic (closer to the satellite image information).

However, an increase in vertical resolution naturally comes with the need for increased computing power, especially when the lowest model level is placed nearer to the surface (as this requires smaller model time steps to maintain numerical stability). Other reasons why the lowest model level cannot be too low are connected with atmosphere/surface coupling. The suggested $\mathrm{MC}\_90$ vertical grid has its lowest model level at around 10 m above the surface (similar to the 65-level grid, which has its lowest level at 12 m) and thus runs with the same model time step as the 65-level vertical grid (75 s). Therefore, the $\mathrm{MC}\_90$ vertical grid requires lower computing power compared to the $\mathrm{MF}\_90$ vertical grid. UWC-West already uses the $\mathrm{MF}\_90$ vertical grid operationally, while MetCoOp is using the $\mathrm{MC}\_90$ vertical grid in pre-operational mode.

5.2. Radiation

The ECMWF Cy25r [14] radiation scheme used in HARMONIE-AROME is no longer maintained and is not up to date with recent developments in radiation schemes currently used in weather and climate models. Moreover, it lacks the flexibility of adding new capabilities, as described in Grailet et al. [96]. In the future, we will use the new ecRad radiation scheme in HARMONIE-AROME. This is the radiation scheme developed at ECMWF and available in HARMONIE-AROME due to sharing parts of the code with IFS. This scheme has been used in the IFS since Cycle 43R3 [97]. Examples of improvements in the scheme are a more recently developed gas optics scheme with updated spectroscopy and the advanced TripleClouds [98] and SPARTACUS [99,100,101] radiative transfer solvers, with the latter accounting for sub-grid cloud 3D radiative effects. The TripleClouds [98] and SPARTACUS solvers were recently optimized and, when combined with the ecCKD [102] spectrally reduced gas optics scheme, were found to be 13× and 2.5× faster, respectively, than the operational IFS radiation scheme [103]. ecRad was implemented in AROME Cycle 49 (another ACCORD CSC), and preliminary tests were carried out.

However, to fully benefit from an improved representation of sub-grid cloud heterogeneity in the radiation scheme, it should probably be run at a coarser resolution. Already at the current horizontal grid step of 2.5 km, the radiative exchange between columns becomes potentially important, and the independent column approximation (ICA), where radiation penetrates only in the vertical, separately within each column, is questionable. For example, Ref. [104] found that local errors in surface SW irradiance associated with ICA become substantial at horizontal resolutions below a few kilometers. To address these issues, a coarse radiation grid (already used in the IFS) is planned for HARMONIE-AROME that will substantially reduce the computational cost while also improving physical realism: information related to sub-grid cloud variability will be aggregated from the relatively fine grid to the coarser radiation grid.

Longer term, developments in fully-3D (inter-column) radiative transfer codes will also be monitored in the hope that their computational cost becomes feasible for high-resolution NWP.

5.3. Cloud Microphysics

In the future, applying the liquid ice multiple aerosol (LIMA) microphysical scheme [105] is planned in HARMONIE-AROME. LIMA, the two-moment bulk scheme, is used in other ACCORD CSCs. It simulates much the same processes and predicts the same hydrometeors as the ICE3 scheme but also includes a prognostic representation of the aerosol distribution. The prognostic concentration of active aerosols is used to calculate the cloud condensation nuclei number concentration, CDNC, ice-freezing and ice crystal nuclei concentrations, and also the concentration of coated ice-freezing nuclei, which are cloud condensation nuclei that have frozen. These variables are, therefore, diagnostic, varying in time, which is an improvement compared with ICE3. Furthermore, LIMA takes more microphysical processes into account than ICE3, such as the Hallett–Mossop process. LIMA can be an alternative option to the combination of NRT aerosols with ICE3 and an option for microphysics in general.

5.4. Scale Aware Shallow Convection

Operational NWP models, with typical resolutions around 2 km, already operate in the gray zone of convection as they (partly) resolve deep convection. With the continuous increase in resolution, more shallow convection will become resolved, and consequently, the transport produced by the convection scheme should be reduced accordingly. Currently, most NWP models deal with the convection gray zone in a very simple way: the shallow convection scheme is shut down if the diagnosed cloud layer depth exceeds a certain threshold. However, even below this threshold, the model can sometimes resolve some convection. In these situations, the convection parametrization consumes too much instability and thereby prevents the model from building up convection by itself. We therefore need more sophisticated ways to adapt the convection scheme to the gray zone.

An approach to make the convection scheme scale aware is presented in [106]. The ideas for this approach originated from the work of [107]. Based on the coarse-graining of LES for several convective cases, the parametrized convective transport is reduced as a function of a non-dimensional grid size. The latter is determined by dividing the grid size by the boundary layer height, where the boundary layer height is used as an estimate of the dominant scale of the convective transport. The larger the scale of the convective transport, the lower the resolution needed to resolve it. In HARMONIE-AROME, an option similar to [106] is available.

As shown by [108], with the boundary layer height for the non-dimensional grid size calculation, the scale of the transport is underestimated in the common case of organized convection. As a result, the reduction in the parametrized convection according to [106], during organized convective conditions, will be too small and will hamper the model’s ability to build up convection by itself. Therefore, we combine the scale aware convection scheme described above with another approach, first presented by [109]. In this approach, the parametrized convection is shut down if the absolute value of the resolved vertical velocity exceeds a certain value, indicating that the model has started to resolve (part of the) convection by itself. Preliminary results in HARMONIE-AROME, using the combined approach of [106,109], are promising.

It should be noted that the mass flux assumptions become unjustified at higher resolutions, and the intrinsically stochastical nature of convection becomes increasingly important. To address the stochastic nature of convection in a simple pragmatic way, a stochastically perturbed parametrization ensemble prediction system can be applied [110].

5.5. Surface

In the future, we plan to switch from the force-restore approach in the ISBA scheme to the multi-layer version of the soil (ISBA-DIF) and the multi-layer snow scheme (ES scheme, the explicit snow scheme) [54]. We expect a better representation of cooling at night with the multi-layer approach, as was shown in preliminary results by ECMWF and Consortium for Small-Scale Modeling (COSMO) [personal communication]. We also hope for a better representation of the annual cycle of soil moisture, especially during spring. We plan to apply the multi-energy-balance (MEB) scheme [111,112], which considers the explicit energy budget for vegetation. MEB contains a detailed description of solar radiation transfer and snow interception in the case of vegetation cover. These schemes are currently under extensive meteorological testing over different domains.

A better representation of turbulent fluxes and wind profiles in the surface layer over high vegetation is expected due to the implemented roughness-sublayer (RSL) scheme [113,114] by Shapkalijevski et al. [115] (Figure 21). This scheme considers the effects of the canopy-induced turbulent mixing in the surface–atmosphere flux–gradient relationships, as well as the dependency of the displacement height and the roughness length on thermal stability over high vegetation. Thereby, it decreases the positive wind bias over these patches (Figure 21b).

Pioneering work has started on producing a new land-cover map with very fine resolution (60 m) over Europe using various sources of information as described in Bessardon et al. [116], Rieutord et al. [117], Walsh et al. [118], Keany et al. [119]. Bessardon et al. [116] used an agreement-based method to combine information from many maps with different semantic, geographical coverage, resolutions, formats, and accuracy to produce a new map called ECOSG+. In addition to ECOSG+, an agreement score map was generated, which can be interpreted as a quality or uncertainty map. Rieutord et al. [117] then used a convolutional neural network to improve the ECOSG+ map in areas where the agreement score map showed low quality and produced a new map called ECOSG-ML. In addition, they suggested producing an ensemble of potential land cover maps with the neural network to reflect the varying quality of the map. Figure 22 shows the ECOSG, ECOSG+, and ECOSG-ML land cover maps for two locations. Further work is needed to fully test these in HARMONIE-AROME, including work on the supplementary physiography datasets of leaf area index, albedo, tree height, and lake depths.

6. General Discussion and Conclusions

In this manuscript, an overview of developments of the HARMONIE-AROME canonical configuration of the ACCORD NWP system between Cycles 40 and 46 is presented. Information on innovations that happened during this lengthy period of time may be of interest to NWP developers from other consortia, forecasters, and academic researchers. Tens of people participated in the model developments. Contributions varied a lot in size and in impact, from small implementations (or even bug corrections, yet having a meteorological impact) to the application of new schemes. Due to that, it is impossible to provide all of the details. We only showed an outline of improvements with some essential features and illustrations and refer the reader to publications about particular developments.

In such a big system as an NWP model, the motivation for particular studies and implementations can be very different. Some developments aim to address well-known model issues like the integral revision of the convection, turbulence, and cloud scheme that improved several aspects of the model but was especially aimed at providing a better representation of low clouds (Section 3.3). Other examples to address important model deficiencies are the calculations of the LW emissivity in Section 3.1 and attempts to improve simulations in stable boundary layer conditions described in Section 4.4. Some developments were devoted to the representation of physical processes, which were missing in the model and believed to be important. In this overview, examples of such studies are the use of NRT aerosols (Section 3.1 and Section 3.2.2) and the implementation of the lake model, FLake (Section 4.7). There are also developments that deal with the implementation of recent knowledge, such as improved physically based methods and algorithms or new databases. Examples of this are the implementation of the wind farm parametrization (Section 3.3.4) and the new physiography land cover map ECOSG (Section 4.1). Moreover, there were efforts made to improve the consistency between different parts of the model, as well as tuning work. The radiation and microphysics schemes were made more consistent through the use of a common cloud droplet number concentration. An example of tuning is the option to include additional 10 m trees to increase the roughness length over low vegetation while using the ECOSG land cover map (Section 4.1).

In the area of model dynamics (Section 2), the most important research and developments described are related to increasing the model resolution, both horizontal and vertical, and improving computational efficiency by using a combination of single and double precision. System aspects usually do not receive much attention in scientific publications. In the area of NWP, however, we think that these are important. A common working environment, facilitated by the model infrastructure, enables the effective exchange of knowledge in a large community and results in an easier application of the code to operational weather forecasting. This infrastructure is built and maintained by the HIRLAM community. A brief description is provided in Section 2.

A significant contribution during the past few years was the development of a method for the coupling and use of external prognostic aerosol data in a limited-area high-resolution NWP model, in particular in the microphysics parametrizations. The methods are described in the sections devoted to radiation (Section 3.1) and cloud microphysics (Section 3.2), while the data sources are described in Appendix A. The methods include: (i) the introduction of aerosol information from external sources (CAMS) in near real-time, (ii) the transport of aerosols by advection, (iii) the description of aerosol effects on radiation processes, and (iv) the description of effects of aerosol on cloud microphysics. We have seen a sensitivity of the resulting radiation and precipitation fluxes to our modifications. However, before operational application, further studies are required to understand the impact of aerosol-related processes and interactions in more detail.

A new, spectral approach for describing LW cloud optical properties was implemented by K.P. Nielsen, (Section 3.1). A comprehensive integral revision of three parametrization schemes, namely, convection, turbulence, and the statistical cloud scheme, described in Section 3.3, resulted in a better representation of especially low clouds. The ICE-T scheme, described in Section 3.2.4, is another outstanding development. It originated from research outside of NWP, namely, from the study of ice loads on transmission lines and aircraft, and resulted in improvements in the representation of supercooled water and the partition between different types of solid precipitation. In the area of surface physics, the most important developments are the implementation of the FLake lake model and a new fine resolution ECOSG land-cover map (Section 4.1).

Prior to accepting new implementations, they are evaluated carefully by developers in specially designed experiments, sometimes using additional observations. The results of some of these evaluations are shown in the paper. The next step is to perform a general evaluation of the model release, which includes all implementations that were added during a certain time period. It is very difficult to show the evolution of general model scores for a long time period covering many implementations. The model is tested for each release for different seasons and domains. In general, for HARMONIE-AROME Cycle 46, the main improvements in the model scores can be seen in the representation of cloudiness. Examples are shown in Figure 23, which demonstrates the improved representation of low cloud cover in HARMONIE-AROME Cycle 46 compared with Cycle 40 over a domain covering Ireland and the United Kingdom during summer and winter. There are improvements in the biases calculated from observations over 140 stations.

Major planned developments in HARMONIE-AROME include evaluating and refining the use of external aerosol concentrations, the ICE-T scheme, and the inclusion of the ecRad radiation scheme and multi-layer ISBA-DIF soil scheme. Work is also on-going on the visibility parametrization, used by the model, as modeling low visibility events is still a challenge. The problem of correctly representing cold temperatures during stable conditions is also a challenge that will need attention in the future.

The HARMONIE-AROME code is not open source, as the ECMWF/IFS code is, in general, not open source. It is open only to the ACCORD community. However, HARMONIE-AROME contains some open source code; for example, the surface modeling platform SURFEX. We can also provide access to the parts of the code that were fully developed by the HIRLAM consortium.