The Cycle 46 Configuration of the HARMONIE-AROME Forecast Model

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 HIRLAM (High Resolution Limited Area Model) international research programme and the ALADIN (Aire Limitée Adaptation Dynamique Développement International) cooperation are introduced. In 2014 HIRLAM and ALADIN agreed on the formation of a single, united consortium, which came to fruition on November 27th 2020. The new consortium, a collaboration between HIRLAM, ALADIN and RC-LACE (Regional Co-operation for Limited Area modelling in Central Europe) is a unique cooperation between 26 countries in Europe and Northern Africa called ACCORD (A Consortium for Convection-scale modelling Research and Development).

There are three model canonical configurations within the ACCORD NWP system. The configuration used by the HIRLAM countries is called HARMONIE-AROME (HIRLAM–ALADIN Research on Mesoscale Operational NWP in Euromed) as introduced in Bengtsson et al. [1]. The so-called HARMONIE scripting system is used to run the AROME (Applications of Research to Operations at Mesoscale) canonical model configuration operationally within the HIRLAM countries. However, several adaptations and improvements have been made to the physical parametrizations in AROME by scientists in the HIRLAM countries - the model configuration which uses these updates is called HARMONIE-AROME. The scripting system facilitates the treatment of the initial and boundary conditions, climate generation, the management and quality control of the observations that are assimilated in the variational analysis process, the execution of the model, postprocessing and the extraction of data for model verification. As well as the scripting system, the HARMONIE software infrastructure consists of a work-flow manager called ecFlow, a work-flow management tool developed by ECMWF (the European Center for Medium Range Weather Forecasting). This software infrastructure is used for both research and operations. The maintenance of such a common working environment has facilitated effective exchange of knowledge resulting in large improvements in weather forecasting.

A single-column model framework, known as MUSC (Modèle Unifié, Simple Colonne) [2], is also available for testing and validating new model developments, and has been used as a validation tool for some of the developments described in this document.

The aim of this technical note is to describe the differences in the reference version of cycle 46 (CY46) of HARMONIE-AROME compared to the cycle 40 (CY40) release described in Bengtsson et al. [1]. We have limited the discussion to the forecast model i.e., the upper-air and surface physics and the dynamics. As well as the reference settings, some non-default research options are also described, components which will be used by default in future cycles of the model. The work is outlined as follows: Changes in upper-air physics since CY40 are detailed in Section 2, surface physics is handled in Section 3 and dynamics and model configuration in Section 4. Upcoming physics developments (upper-air and surface) are discussed in Section 5.

2. Upper-Air Physics

2.1. Radiation

The Morcrette radiation scheme from ECMWF’s IFS (Integrated Forecasting System) cycle 25R1 is the default shortwave (SW) radiation scheme in HARMONIE–AROME and contains six spectral intervals. The default longwave (LW) radiation scheme uses the Rapid Radiative Transfer Model (RRTM) of Mlawer et al. [3] and contains 16 spectral bands. Both the SW and LW schemes are described in the IFS [4] and the mesoscale research model Meso-NH [5] documentation. Because of computational constraints the full radiation calculations are currently performed every 15 minutes. An alternative broadband radiation scheme called ACRANEB2 [6,7] is also available in the model, though not used operationally by default. Comparisons over Finland have demonstrated a general agreement of the clear-sky SW and LW radiation fluxes at the surface between the default radiation scheme and ACRANEB2 [8,9]. Further details on the use of these schemes in HARMONIE-AROME are described in Bengtsson et al. [1] with new developments since then described here, in particular the following three changes: 1) the inclusion of Copernicus Atmosphere Monitoring Service (CAMS) near-real time (NRT) aerosols, 2) changes regarding the cloud droplet number concentration, and 3) changes to the LW cloud liquid optical property parametrization.

The radiative transfer calculations use the following inherent optical properties (IOPs) of cloud particles, aerosols and atmospheric gases: optical thickness, single scattering albedo (SSA) and asymmetry factor (g). Regarding cloud particles we mean prognostic specific cloud liquid and cloud ice content. The aerosols consist of monthly climatologies and the atmospheric gases consist of prognostic H${\hspace{0.17em}}_2$O, a fixed composition mixture of CO${\hspace{0.17em}}_2$, N${\hspace{0.17em}}_2$O, CH${\hspace{0.17em}}_4$, and O${\hspace{0.17em}}_2$ and monthly climatologies of O${\hspace{0.17em}}_3$). By default, the Tegen aerosol climatology [10]) consisting of vertically integrated aerosol optical depths (AOD) of land, sea, desert, and urban tropospheric aerosols at a wavelength of 550 nm (AOD550) are used in the radiative transfer calculations, along with prescribed constant background tropospheric and stratospheric aerosols.

The first new radiation-related development described here is the use of NRT aerosol concentration data, provided by CAMS as suggested by [11]. These data are now included in both for radiation and cloud microphysics parametrizations (Section 2.5), while in the default system the climatological AOD550s only influence the radiative transfer meaning that the use of aerosols in radiation and microphysics is not consistent in the older configuration. Three-dimensional NRT aerosol mass mixing ratio (MMR, unit kg kg${\hspace{0.17em}}^{-1}$) fields are introduced via the first guess and lateral boundary conditions (LBC) of the model, and advected by the model dynamics. 14 aerosol species are included: 3 sea salt species (fine, jet and spume drop modes), 3 desert dust (fine, coarse and supercoarse) 2 organic matter (hydrophilic and hydrophobic), 2 black carbon (hydrophilic and hydrophobic), one Sulfate, 2 Nitrate (fine and coarse mode) and Ammonium. All aerosols can be removed by wet deposition, while the coarser modes can also be removed by dry sedimentation. With the exception of sea salt emissions, aerosol sources or generation mechanisms are not yet implemented in HARMONIE-AROME.

Three-dimensional AOD550 fields are obtained for the radiation calculations from the NRT aerosol MMR fields using the mass extinction (ME) coefficients (provided by [12]) for 11 species. A relative humidity of 80 % is assumed for conversion of the MMRs to AOD fields, AOD550(x,y,z,species) = MMR(x,y,z,species) x ME(humidity,wavelength,species). 11 species are used, as opposed to the 14 mentioned in the previous paragraph, because currently the nitrate and ammonium aerosols are not included in the radiation calculations. The AOD550s are grouped into the same general categories as the Tegen species, to enable these to be readily used by the default radiation scheme. In addition to ME, the values of the single-scattering albedo and asymmetry factor and the spectral distributions of all aerosol inherent optical properties (AIOPs) are, as previously, prescribed. Note also that the cloud optical property parametrizations [1] remain unchanged when aerosol MMRs are introduced.

The impact of aerosols on SW radiation was found to be large in the case of desert dust intrusions as shown in Figure 1 and Figure 2. Generally in clear sky cases, SW radiation was found to be slightly underestimated with NRT aerosols compared to when the Tegen climatology is used. Figure 1 shows a mostly clear-sky dust example that occurred on the 20th of February 2023 over the Iberian Peninsula, where Saharan dust reduced the global SW radiation by nearly 100 W/m${\hspace{0.17em}}^2$ at some stations. The daily mean global radiation for Tegen (REFERENCE) and CAMS NRT aerosol experiments is shown in Figure 1. Hourly global radiation is shown in Figure 2 compared with observations, where the global radiation forecast is overestimated by the REFERENCE Tegen experiment, while the introduction of NRT aerosols reduces the global radiation to values similar to those observed.

The second new radiation-related development described in this section concerns cloud droplet number concentration. Radiation schemes are known to be sensitive to the definition of cloud ice and water particle sizes. The optical density of a cloud of a given total water content distributed to a multitude of small droplets is significantly larger than that of a cloud consisting of fewer larger droplets. The cloud droplet number concentration (CDNC) can be used to calculate the cloud liquid droplet effective radius (r${\hspace{0.17em}}_{\mathrm{eff}}$) for the radiation scheme. In CY46, for consistency reasons, the CDNC used in the radiation and microphysical schemes is the same; this wasn’t the case in CY40, where the prescribed continental and marine aerosol coefficients were different in both schemes. A vertical profile of CDNC is used now in place of the coefficients and is described in Section 2.5 along with other details.

The remainder of this section is dedicated to the third radiation-related development concerning LW cloud liquid optical properties. Long-standing issues forecasting fog in HARMONIE-AROME (fog was too dense, too cold and too widespread) led to a Masters thesis by Tosca Kettler in collaboration with KNMI [13]. Her study, primarily using MUSC, investigated the strong cooling of the fog layer and a possible link with LW optical depth. By default, the [14] LW effective emissivity, $\epsilon$, was used in HARMONIE-AROME as described in Equation 1:

$$\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), 0.144 is an empirically fitted coefficient, and 1.2 and -0.006 are coefficients that describe the linear dependence on the effective radius ($r_e$). 1.66 is the diffusivity factor [15]. The 0.144 coefficient 1 was shown to be too high when [13] used data from the Cabauw supersite to derive a more realistic LW effective emissivity as shown in Figure 3 (coefficient of -0.096 instead of -0.144). This effective LW emissivity was derived from downward LW radiative fluxes at the top of the Cabauw tower and at the surface. An empirical relationship was then derived between the effective emissivity and the LWP calculated from visibility measurements at 6 levels in the Cabauw tower. In addition to the LW radiation, the cloud droplet number concentration was found to be very important regarding the fog problems (see Section 2.5 for more details).

Theoretical analyses using Mie computations for the LW spectrum confirmed that the HARMONIE-AROME LW cloud emissivity was indeed too high. The results of these calculations for the 16 LW spectral bands in the radiation scheme are shown in Figure 4. Here, it can be seen that the Kettler parametrization fits the theoretical results well for wavelength bands 8, 9 and 10 and the smallest, and most common, cloud droplet sizes. These spectral bands are the most important for the cloud LW radiative effects, as they are least affected by greenhouse gas absorption. However, in an atmosphere with highly variable water vapour concentrations, the important wavelength bands vary. Thus, a LW cloud emissivity parametrization using a single coefficient is not ideal. A new LW cloud liquid optical property scheme, the Nielsen scheme (personal communication), was developed, where the 16 LW spectral band cloud liquid droplet mass absorption coefficients were parametrized based on Mie-Debye theory. This parametrization is given in Equation 2, where n(r${\hspace{0.17em}}_e$) is the spectral band mass absorption coefficient, r${\hspace{0.17em}}_e$ is the cloud droplet effective radius and a${\hspace{0.17em}}_{\lambda }$-g${\hspace{0.17em}}_{\lambda }$ are the coefficients for each LW spectral band, as detailed in Table 1.

$$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\text{µ}\mathrm{m}$$

(2)

2.2. Convection

[16] did a comprehensive integral revision of three parameterisation schemes in the HARMONIE-AROME model that collectively parameterise boundary layer processes: the cloud scheme, the turbulence scheme, and the shallow cumulus convection scheme. These updates mean that the CY46 cloud, turbulence and convection schemes differ significantly from the implementation described in [1]. An insight to these model updates on subgrid-scale responses, as well as on the grid-scale, for a cold air outbreak (CAO) case is described in [17]. In this section we focus on convection, where the biggest change involves the coupling of the convection and the turbulence schemes.

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 parameterize deep convection. However, shallow convection still needs to be parameterised. The shallow convection scheme in HARMONIE-AROME applies a dual mass flux framework [18] in which two updrafts are distinguished: a dry updraft that does not enter the cloud layer as a saturated updraft, and a moist updraft that reaches the lifting condensation level (LCL) and continues its ascent in the cloud layer. The shallow convection scheme is described in detail by [16] where all the modifications compared to the CY40 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 [19] is used as a source term in the turbulent kinetic energy (TKE) equation. This particularly enhances the subcloud-to-cloud layer transport in better correspondence with Large Eddy Simulation (LES) results for shallow convection - see [16] and Figure 5 which shows the total turbulent transport of moisture $\left({\overline{w^{\prime }{r}^{\prime }}}_{t,tot}\right)$ according to the DALES LES model, HARMONIE-AROME with an energy cascade and HARMONIE-AROME without an energy cascade. In the case of HARMONIE-AROME, ${\overline{w^{\prime }{r}^{\prime }}}_{t,tot}$ is the sum of the transport by the turbulence and convection schemes, whereas for the LES 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 5 reveals that like the LES model, the version with the energy cascade term describes the increasing moisture transport with height (i.e., drying) in the sub-cloud layer. This results in a much better representation of the ventilation of the sub-cloud layer and moistening of the cloud layer.

2.3. Turbulence

In this section we focus on the turbulence scheme where changes have been made to improve the forecasts of low cloud. Another turbulence-related update since CY40 is the inclusion of a wind farm parametrization in the model.

HARMONIE-AROME releases since CY36 use the HARATU (HArmonie with RAcmo TUrbulence) [1,16,20] turbulence scheme in place of the CBR (Cuxart-Bougeault-Redelsperger) scheme [21]. Both schemes combine a prognostic equation for turbulent kinetic energy (TKE) with a diagnostic length scale. Although the HARMONIE-AROME model significantly improved on several aspects of the turbulence scheme with the transition to HARATU (see [1,22,23]), HARATU also contributed to an underestimation of low cloud cover and an overestimation of cloud base heights. In the comprehensive integral approach mentioned earlier the turbulence scheme was revised to substantially improve on forecasts of low clouds in particular (see Figure 6). This figure shows the frequency bias of cloud base height classes for one month using the CY40 reference [1] (left panel), and with the modifications of [16] (right panel). The blue, green, and orange lines refer to +3, +24, and +48h forecasts, respectively. In the reference CY40 version $+24h$ forecasts contain fewer than 20% (y-axis) of the observed number of cloud base heights of around 178 feet (x-axis). The updated version on the right panel clearly improves the climatology of low cloud base heights. Two modifications to the turbulence scheme are most relevant to the improved low cloud climatology. The first involved a re-tuning of the scheme based on the Monin–Obukhov similarity theory, following [24] and [25]. The second involved 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 but, for example, the triggering of deep precipitating convection is also influenced [16].

Figure 6. Frequency bias of the cloud base height in feet (1 ft is 0.3048 m) for December 2018 with (a) CY40 [1] and (b) CY40 with the updates in [16] as applied in CY43 and CY46. The blue, green, and orange lines refer to +3, +24, and +48h forecasts, respectively. Note the extreme underestimation of cloud bases around 178 ft (approximately 54 m) in CY40. Fewer than 20% of the observed number of cases are actually predicted in +24h CY40 forecasts, with a clear improvement seen when the latest changes were introduced. European Geosciences Union 2022, from Figure 20 of [16].

Figure 6. Frequency bias of the cloud base height in feet (1 ft is 0.3048 m) for December 2018 with (a) CY40 [1] and (b) CY40 with the updates in [16] as applied in CY43 and CY46. The blue, green, and orange lines refer to +3, +24, and +48h forecasts, respectively. Note the extreme underestimation of cloud bases around 178 ft (approximately 54 m) in CY40. Fewer than 20% of the observed number of cases are actually predicted in +24h CY40 forecasts, with a clear improvement seen when the latest changes were introduced. European Geosciences Union 2022, from Figure 20 of [16].

The second important development regarding turbulence involves the introduction of a wind farm parametrization (WFP). As investment in renewable energy is increasing, wind turbines are increasingly a part of the European landscape. In order to model the effect that these (large) wind turbines have on the atmosphere a WFP was introduced in HARMONIE-AROME. Van Stratum et al. [26] described the initial implementation of the WFP and showed an evaluation using one year of HARMONIE-AROME simulations. The parameterisation 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), where U is the two wind components (u and v), $|\overrightarrow{V}|=\sqrt{u^2+v^2}$, A is the rotor area, and ${\Delta }_k$ is the volume of the grid cell. The WFP requires thrust coefficients ($C_T$) and power coefficients ($C_P$) for each wind speed and each wind turbine. The remainder of the energy extracted from momentum, that is not converted to power, is transferred to turbulent kinetic energy, with $C_{\mathrm{TKE}}$ = $C_{\mathrm{T}}$ - $C_{\mathrm{P}}$.

2.4. Statistical Cloud Scheme

In HARMONIE-AROME, ice clouds, which are mainly high clouds, are parameterized separately in a relative humidity scheme (see Section 2.5.3). On the other hand, water and mixed-phase clouds, which commonly appear as low and medium clouds, are estimated with a statistical cloud scheme in which the cloud fraction and liquid water content can be derived once the sub-grid variability of temperature and moisture is known [27]. A detailed description of the statistical cloud scheme, including several fixes and modifications to the CY40 version, can be found in [16]. The main modifications concern:

• A proper derivation of the thermodynamic coefficients.
• The inclusion of the covariance term of temperature and humidity.
• Using a dissipation length scale that is consistent with the one in the turbulence scheme.
• The description of the dissipation of the variances is more consistent with literature.
• An erroneous factor 2 for the convective contribution to the variance is removed.

Although the updated statistical cloud scheme is based on a more sound physical basis and several errors are removed, results for the idealised ARM cumulus case [28] revealed an even more pronounced underestimation of the variance, see Figure 7. The variance is especially underestimated in the upper part of the cloud layer, leading to an underestimation of the small cloud fractions at these heights [16]. To solve this problem, we are currently investigating a physically plausible approach in which the current constant dissipation time scale of convection in the cloud scheme is made height dependent.

2.5. Cloud Microphysics

The microphysical core in HARMONIE-AROME is a one-moment bulk scheme containing three different classes of ice parameterization based on developments originally done in Meso-NH [29,30]. This classification is commonly referred to as "ICE3". These three classes of ice, the solid hydrometeors, are snow, cloud ice and a combination of hail and graupel. In addition, three liquid or gaseous hydrometeors of rain, cloud water and water vapour are taken into account. All of the these hydrometeors are represented by prognostic mixing ratios and advected by the model dynamics; horizontally by the semi-Lagrangian scheme [31] and vertically by the sedimentation process [32]. The grid box-wise particle sizes are estimated from a generalized Gamma distribution.

This section on cloud microphysics is split into four subsections covering the following developments: 1) changes to the CDNC important for fog and clouds in general, 2) the use of NRT aerosols in microphysics, 3) the OCDN2 option in the HARMONIE-AROME version of ICE3, and 4) ICE-T which focuses on improving the representation of supercooled liquid (SCL).

2.5.1. Cloud Droplet Number Concentration

As described in more depth in Section 2.1, by default HARMONIE-AROME does not predict aerosol concentrations in the forecast. Therefore, microphysically important variables such as the number concentration of cloud condensation nuclei, cloud droplets, and ice crystals, are prescribed by default. The prescribed values are either height dependent, as described in Contreras Osorio et al. [33] or approximated process-wise.

By default, the CDNC in HARMONIE-AROME CY46 now has a vertical dependence on height, with the same profile used over land and sea areas. The reason for using the same profile was to eliminate the artificial reduction in stratiform precipitation that used occur at the land sea boundaries due to the use of different CDNCs. A reference concentration of $250 c{m}^{-3}$ is considered at a pressure of 1000 $hPa$. The concentration increases linearly with pressure at a rate of $P/1000$, where P is the pressure at each model level. A reduction to CDNC is applied at the lowest model level, where the concentration is multiplied by a factor of $0.25$ - this was done to improve the forecasting of low visibility and fog because using the same CDCN profile over land and sea led to too much fog over the sea areas. Prior to CY43 constant values of CDNC were used for all model levels, with a value of 500 cm${\hspace{0.17em}}^{-3}$ for urban areas, 300 cm${\hspace{0.17em}}^{-3}$ for land areas and 100 cm${\hspace{0.17em}}^{-3}$ over the sea/ocean. These values were shown to result in an over-estimation of fog, and an overestimation of cloud condensate in the lowest thickest clouds.

An example of the issue with the condensate in low thick clouds is shown in Figure 8 where an MSG visible satellite image for 12 Z on July 8th 2019 is shown along with the MSG Seviri cloud water path product from KNMI and output from two HARMONIE-AROME experiments, one with default CY43 settings (CDNCs of 500 cm${\hspace{0.17em}}^{-3}$, 300 cm${\hspace{0.17em}}^{-3}$ and 100 cm${\hspace{0.17em}}^{-3}$ over urban, land and sea areas respectively) and the other with a CDNC of 50 cm${\hspace{0.17em}}^{-3}$ everywhere and at all vertical levels (this was later replaced by the profile discussed above), and the LW effective emissivity suggested by Kettler [13]. It is clear from Figure 8 that the cloud water path was overestimated in the default CY43 using the old CDNC and LW effective emissivity. A CDNC of 50 cm${\hspace{0.17em}}^{-3}$, with a re-tuned LW effective emissivity gave much better results. The most recent configuration in CY46 (profile of CDNC and the Nielsen 2020 LW cloud liquid optical property scheme) also shows significant improvements in the cloud water path (not shown). Further experiments led to the use of a CDNC profile as described above and the Kettler [13] LW effective emissivity tuning was replaced by the more robust Nielsen LW cloud liquid optical property scheme.

2.5.2. Use of CAMS NRT Aerosols in ICE3

A configuration to use the three-dimensional aerosol MMRs in the cloud-precipitation microphysics was introduced by [11] but is not yet a default option. See Section 2.1 for the description of the aerosol species that are introduced in the forecast model for both radiation and cloud microphysics parametrizations. The key variable derived from the aerosol fields is the CDNC. Estimations of CDNC are 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 thermodynamical variables and the vertical velocity. Furthermore, the CDNC is used in the parametrization of various processes leading to the growth of activated cloud droplets to liquid and solid precipitation, as discussed in Section 2.5.3 and Section 2.5.4.

Case studies [11,34] show that precipitation, in particular the phase of the precipitation - snow, graupel, rain, is affected by the introduction of aerosol MMRs in the cloud microphysics. Changes in cloud distribution and optical thickness lead to changes in the surface radiation fluxes in cloudy cases.

Figure 9 shows 1D histograms of observed versus modelled clear sky index (CSI) for Tegen and CAMS NRT aerosol experiments, where the clear sky index is the ratio of the clear sky global SW radiation and the observed or modelled global SW radiation (the lowest values thus correspond to cloudy conditions, with values close to 1 representing clearer skies). SW radiation observations from 20 synoptic sites around Ireland were used, along with corresponding model data for the same locations. The results for a 2 week Summer period (June 1st to 14th 2018) are shown on the left, with a Winter period (February 3rd to 17th 2020) shown on the right. A clear overestimation of low CSI in the Tegen experiment can be seen in both seasons, consistent with an overestimation of cloud condensate in the thickest clouds (shown in Figure 8). This overestimation of low CSI is not present when CAMS NRT aerosols are used.

In the case of CAMS NRT aerosols a general overestimation is seen in global SW radiation in Figure 10, where the positive and negative biases are plotted on the positive axis to highlight whether an experiment results in more positive or negative biases overall. These biases were calculated for the same 20 stations as before. It is clear that in the Tegen experiments (both Summer and Winter; left figures) there are negative biases overall, with positive biases overall when CAMS NRT aerosols are used (right figures). Note that the differences arise from the impact of the aerosols in both the radiation and the microphysics schemes.

2.5.3. OCND2

As described in [1] and [35] ICE3 was supplemented with an option called OCND2 due to weaknesses found in relation to stable boundary layer conditions over Northern Europe, where the model generated ice too quickly when supercooled liquid was expected. Spurious fog was also produced at temperatures below -20${\hspace{0.17em}}^{\circ }$C. The major reason for such systematic model deficiencies is linked to the treatment of mixed-phase and pure ice clouds. As well as with the OCND2 option, the treatment of supercooled liquid was also improved by upgrading ICE3 to ICE-T as described in SubSection 2.5.4. In particular, the following aspects of the ICE3 microphysics parametrization were modified under OCND2 in order to address weaknesses:

• The separation between liquid water processes and ice water processes was improved. This means that the statistical cloud scheme (See Section 2.4) only deals with cloud liquid water, including cases when temperatures are below freezing. Thus, all ice processes are taken care of by the OCND2 version of the ICE3 scheme.
• Evaporation/deposition of cloud ice water is a conversion between ice and vapour and not between ice and liquid.
• The deposition rate of the ice water species was reduced.
• The cloud cover, from the point of view of users of the forecast (the public), was modified to account for the lower optical thickness of ice clouds compared to water clouds.
• The ice number concentration was reduced between temperatures of $0{ }^{\circ }$C and $-25{ }^{\circ }$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
• Avoid calculations of saturation pressure when the saturation pressure is near or above atmospheric pressure. This is done just for technical reasons, and affects calculations in the Stratosphere only.
• In order to save computing time, the ICE3 scheme should be active only when any non-vapour water species are present above a low threshold, or when the air temperature is below freezing. Unfortunately, this did not always happen when the second criteria was satisfied. This occasionally led to spurious square-like ice-clouds where areas with sufficient water species are surrounded by areas with too little cloud ice water, as shown in Figure 11. A fix has now been implemented.

With OCND2, the microphysics scheme is more "liquid friendly", which means that supercooled droplets do not freeze as quickly. However, it is still not "liquid friendly" enough. There are cases where supercooled clouds seem to disappear too quickly. Therefore, an extension of the OCND2 scheme labeled "ICE-T" has been developed and is described in SubSection 2.5.4.

2.5.4. ICE-T

The development of ICE-T was motivated by the estimation of ice loads on transmission lines from atmospheric icing. When using the supercooled liquid water (SCL) from HARMONIE-AROME simulations for the estimation of ice loads, it was found that, despite the efforts made with "OCND2", the model still had a tendency to glaciate the clouds prematurely, and hence underestimate ice loads on transmission lines [36]. Elements from the Thompson microphysics scheme [37] found in the Weather Research and Forecasting (WRF) model, were implemented in the ICE3 microphysics with the OCND2-option active. The name of the new option is called ICE-T, and reflects the combination of ICE3 and the Thompson scheme. The changes are described in detail in [38], but the most prominent ones 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 [36] and aircraft [39]. The results show a clear shift towards more SCL water. However, the atmospheric content of ice species is also increased. This is due to a shift from graupel to snow, and since snow has a slower terminal fall speed than graupel, the residence time, and hence accumulated amount, is increased. There is an increase in surface precipitation as snow, and a decrease in 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.

Both studies found that the increased SCL was in better correspondence with observations of ice loads, measurements of atmospheric content of liquid and ice water by satellites, and pilot reports of experienced aircraft icing. The shift in the precipitation pattern is less beneficial and needs to be explored further.

ICE-T was tested during a helicopter measurement campaign launched from Alta, Norway, in April 2023. The campaign was lead by Airbus in order to test their helicopters’ ability to fly through heavy icing conditions. On April 19th 2023, heavy icing conditions inside lenticularis clouds occurred over the mountainous areas in the vicinity of Alta. The helicopter was equipped with a Cloud Droplet Probe (CDP), which measures hydrometeors in the range of $2-50\mu \mathrm{m}$, which are essentially cloud droplets. During an afternoon flight they measured liquid water content of mostly 0.8 gm${\hspace{0.17em}}^{-3}$ and above. The highest value measured was about 1.3 gm${\hspace{0.17em}}^{-3}$.

Two parallel simulations were carried out using CY46, one with the default ICE3 and the OCND2 switch active, here called DEF, and one where ICE-T was active, called ICE-T. The simulations were initialised on April 19th 00:00 UTC, with no upper-air data assimilation. Figure 12 shows cloud water content for both simulations and the difference between them, at model level 41 corresponding to approximately 820 hPa, where the helicopter encountered the heaviest icing conditions. Overall, the ICE-T simulation has higher cloud water content than DEF, and is closer to what was observed by the CDP.

Figure 12. Simulated cloud liquid water content in the Alta region for model level 41 (approximately 820 hPa) at 14:00 UTC April 19th 2023 for DEF (left), ICE-T (middle) and the difference (right).

Figure 12. Simulated cloud liquid water content in the Alta region for model level 41 (approximately 820 hPa) at 14:00 UTC April 19th 2023 for DEF (left), ICE-T (middle) and the difference (right).

3. Surface Physics

In HARMONIE-AROME, the SURFEX modeling platform [40] is used to model the surface of the Earth. This platform is developed and maintained mainly by Météo-France. In CY46 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 4 main types of surface (so-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 as 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 facilitates the exchange of turbulent fluxes. The nature tile can be further divided into sub-types (so-called patches) for different kinds of vegetation, e.g., crops, forestry, rocks (with no vegetation). SURFEX also includes the software needed to initialise models and prepare land-cover information, including all of the necessary external parameters such as leaf area index, albedo, tree heights and lake depths. For the nature tile, the vegetation parameters are specific to each patch. Fluxes and diagnostics are calculated for all tiles and patches and then aggregated for each grid box according their corresponding fractions.

The following subsections outline developments done 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 in the sea and inland water tiles.

3.1. Physiography

The land cover characteristics in HARMONIE-AROME CY46 are described using the ECOCLIMAP Second Generation (ECOSG) land cover map (https://opensource.umr-cnrm.fr/projects/ecoclimap-sg) developed by Météo-France, which has a primary resolution of ca. 300 m. This is an advancement compared with CY40, when the 1-km resolution ECOCLIMAPv.2.2 [41] was used.

ECOSG was corrected for use in HARMONIE-AROME CY46 using the Randolph Glacier Inventory (https://www.glims.org/RGI/) for 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. To represent topography, the GMTED2010 dataset (https://topotools.cr.usgs.gov/gmted_viewer/viewer.htm) is used. It is corrected using data from the ArcticDEM (https://www.pgc.umn.edu/data/arcticdem/) and Copernicus GLO-90 (https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.1) projects. For the soil properties, we rely on the SoilGrids database (https://www.isric.org/explore/soilgrids) which has a primary resolution of ca. 300 m. For lake depths, we use the Global Lake Database v.3 [42], with a primary resolution of ca. 1 km, and for the initialization of lake variables during the first model run we use the global lake climatology by [43].

3.2. The Urban and Nature Tiles

The urban tile is simulated by the Town Energy Balance (TEB) model [44]. 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 in SURFEX. Currently, HARMONIE-AROME CY46 uses the force-restore approach (ISBA-FR) [45] version of the scheme. To describe the evolution of snow cover, the one layer snow scheme [46] is applied. Diagnostic temperature, humidity and wind at 2 and 10 m over land are calculated by means of interpolation in the surface layer, with a dependency on turbulent fluxes. In CY46 we increased the number of nature patches from 1 to 2 so that we now have distinct low and high vegetation patches.

While using 2 patches in combination with the high-resolution ECOSG land cover map, we experienced an increase in 10 m wind speeds 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 the 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 it 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 of parameterizing this effect by increasing the vegetation height is included in HARMONIE-AROME CY46 (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. LFAKETREE 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 13 shows verification scores for two test periods, with different physics options: CY40 (red), CY43 default (green), CY43 with "LFAKETREE" (blue). Only results from the 00:00 UTC forecast cycles are shown, in order to highlight any diurnal pattern. From Figure 13 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 ECOSG tree height values led to roughnesses that were too high, and hence an underestimation of surface wind speeds. To circumvent this problem, a tree height scaling with latitude was introduced in HARMONIE-AROME.

3.3. Snow Melt Adjustment

The one-layer snow scheme [46] has a problem with melting very small snow amounts. This issue is caused by the use of 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^{2}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 doesn’t 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 14.

3.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 [47,48]. Too much heat and momentum transport results in a warm bias in near-surface temperatures and reduced low-level wind speeds [49]. Halting this transport triggers a runaway cooling regime that causes unrealistically cold surface temperatures [50,51]. Turbulent flux calculations in the SBL are frequently the target of model tuning to enhance performance. Such tuning efforts have also continued using HARMONIE-AROME. The main tuning parameter is a threshold for the gradient Richardson number. In the model there is an upper limit for this threshold, so that the calculated gradient Richardson number can’t be larger than it. This leads to a limitation in the calculations of the drag coefficients and hence the turbulent fluxes.

The default maximum value of the gradient Richardson number (XRIMAX) in the previous version of HARMONIE-AROME was 0.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, but 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 [52] suggested modifying the dependency of the drag coefficients on the gradient Richardson number, restricting the drag coefficients to be equal to neutral values, for a certain interval of stability characterized by the gradient Richardson number. Model performance is very sensitive to the different options used for calculating the drag coefficients and values of the maximum gradient Richardson number. However it is impossible to find an optimal solution, which would result in an equally good performance for all geographical conditions.

Algorithms used to interpolate meteorological variables to screen level values can differ, which also affects the verification scores. Screen-level 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. For a detailed depiction of the impacts of different options on stable boundary layer behaviour the reader is referred to Kähnert et al. [53]. The default maximum value of the gradient Richardson number in the reference version of HARMONIE-AROME CY46 is 0.2 but this value does not give the best scores in domains dominated by snow in winter.

3.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 [54], possibly with directional components based on Georgelin et al. [55]. 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, and the lowest model layer (currently ca. 14 m). This may lead to nonphysical behaviour 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 the airflow over hills and mountains [56] in the form suggested by Beljaars et al. [57].

A third option called OROTUR was suggested and implemented according to [58]. This is another realization of the Wood et al. [56] idea. OROTUR 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 lowest model level wind speed 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 to $\Delta x^2$ allows the orography 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 no tiles or patches are considered. The method includes two tunable constants. The default version of HARMONIE-AROME CY46 does not use any small-scale orography parametrization. However in the Spanish Meteorological Agency (AEMET) OROTUR is applied operationally. [59] reported that the scheme contributed to a reduction in the positive bias in 10 m mean wind speed and gusts. In HARMONIE-AROME, vertically propagating buoyancy waves from orography are not parametrized because the generation of these waves is believed to be described by the fine-resolution non-hydrostatic dynamics of the model. Dissipation of the waves should be taken care by the turbulence parametrizations, as suggested by Bougeault and Lacarrere [60].

3.6. The Sea Tile

For use in 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 turns uses the OSTIA analysis (https://ghrsst-pp.metoffice.gov.uk/ostia-website/index.html) with some modifications. Other SST sources are also used; for example MetCoOp 2 uses the results of the NEMO ocean model [61] 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, e.g., Belušić et al. [62].

To calculate turbulent fluxes in the surface layer over unfrozen seas, the Exchange Coefficients from Unified Multicampaign Estimates (ECUME) scheme by [63] is used. SURFEX provides two versions of this scheme, namely ECUME and ECUME6 [40]. They were developed using observations from different measurement campaigns and use slightly different parameters and formulations. Extensive testing was performed in order to assess the impact of these schemes. Observations from the EUREC⁴A project (https://eurec4a.eu/) were used in the evaluation process. Comparisons with the EUREC⁴A observations were performed in phase-space, see [64] for explanations. Figure 15 shows some results from this study, where the ECUME6 scheme overestimates the latent heat flux (LHF). However, the bias is small for the ECUME scheme. These results were also confirmed by climate run validation (Oskar Landgren and Bert van Ulft, personal communication). The default scheme in the reference version of HARMONIE-AROME CY46 is therefore ECUME.

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

3.7. The Inland water Tile

For the inland water tile (lakes, reservoirs and rivers) HARMONIE-AROME uses the FLake [67] lake model. FLake is parametric: it uses a self-similarity concept and predicts characteristics of the temperature profile in the water column, ice and snow. This approach incorporates much of the essential physics and offers a good compromise between physical realism and computational cost. The temperature profile in the water column is represented by the mixed layer and the thermocline. The mixed layer depth, the mean water temperature and the bottom temperature are predicted. The integral of the temperature profile in the thermocline (the so-called shape factor) is also predicted. The water surface temperature is diagnostic. The temperature profile in snow and ice are assumed to be linear. The model predicts the depths of snow and ice and the snow and ice temperatures. The snow block uses modifications by [68]. 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 (assuming a certain error for saline lakes).

4. Dynamics and Model Configuration

The dynamical core of HARMONIE-AROME was previously described in Bengtsson et al. [1], and just a brief overview of the default configuration for CY46 is given here.

The nonhydrostatic core was developed by ALADIN [31,69] and is shared with the ALADIN and AROME [70,71] canonical model configurations. It uses the fully compressible Euler equations with a mass-based hybrid pressure terrain-following vertical coordinate [72]. 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 [73].

The time scheme allows for a general iterative centred implicit (ICI) approach for improved stability [31]. Typically ALADIN and AROME use a predictor–corrector method in operational configurations. HARMONIE-AROME, on the other hand, uses a single-step SISL discretisation (note that in practice the two approaches have a comparable cost, as the predictor–corrector allows for a slightly longer time step which off-sets the additional effort). Stability is enhanced through the use of the SETTLS [Stable Extrapolation Two-Time-Level Scheme, [74], 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 [75] for the upper boundary condition has further helped to stabilise integrations.

The default, reference configuration of HARMONIE-AROME uses a horizontal grid-spacing of 2.5 km, with a 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, and 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 Bengtsson et al. [1], there is an option to run 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 operations, with minimal impact on accuracy.

Linear spectral horizontal diffusion is applied to the spectral prognostic variables. In addition, the flow-dependent nonlinear SLHD scheme [Semi-Lagrangian Horizontal Diffusion, [76] is used for hydrometeors and cloud water. Similar to as described in Seity et al. [71], quasi-monotonic interpolation is used in the semi-Lagrangian advection, while the COMAD scheme [Continuous Mapping about Departure points, [77] is also employed for the dynamical fields.

By default, HARMONIE-AROME runs in double precision (i.e., it uses 64 bit variables for storing floating-point data, as opposed to a single precision configuration which uses 32 bit floating-point variables), however since Cycle 43h2 single precision forecasts can be activated by using the “dual” precision option. In this case, the forecast model alone runs in single precision, while all other model components are retained in double precision, which ensures that highly precision-sensitive components, such as data assimilation, are not impacted. This is the approach adopted at ECMWF where single precision forecasts, providing a runtime saving of approximately 40% without negatively impacting forecast quality, have been used operationally since Cycle 47r2 [78].

Limited testing of this option in HARMONIE-AROME CY46 has observed runtime savings of about 35% with little to no impact on point verification scores, apart from a small positive mean sea level pressure bias in single precision relative to double. This bias can grow with lead time and is most likely related to radiation issues (i.e., the air cools too much in single precision). More extensive testing of dual precision has been carried out with Cycle 43, where single precision forecasts are now used operationally in the United Weather Centres – West consortium 3.

5. Upcoming Developments in HARMONIE-AROME

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

5.1. Radiation

The Morcrette scheme used in HARMONIE-AROME is by now severely outdated. Future cycles of HARMONIE-AROME are expected to use ecRad instead, the modular ECMWF radiation scheme used in the IFS since cycle 43R3 [79]. This should represent a major improvement by enabling the use of, for example, modern gas optics schemes with updated spectroscopy and more advanced radiative transfer solvers such as SPARTACUS, which accounts for sub-grid cloud 3D radiative effects [80]. The TripleClouds [81] and SPARTACUS solvers were recently optimized, and when combined with the spectrally reduced gas optics scheme ecCKD [82], found to be 13× and 2.5× faster respectively than the operational IFS radiation [83]. This version of ecRad has been implemented in a CY49 branch of AROME and preliminary tests have been carried out.

However, to fully benefit from an improved treatment of sub-grid cloud heterogeneity in the radiation scheme, it should probably be run at a coarser resolution. Already at the current horizontal grid-spacing of 2.5 km, radiative exchanges between columns become potentially important, and the Independent Column Approximation (ICA, whereby radiation is treated in 1D separately within each column) is questionable. For example, [84] 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, which will substantially reduce computational cost while also improving physical realism: ecRad inputs related to sub-grid cloud variability can then be computed from the fine grid relative to the radiation grid, as opposed to using constant values.

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.2. Scale Aware Shallow Convection

Operational NWP models, with typical resolutions around 2km, already operate in the gray zone of convection as they (partly) resolve deeper convection. With the continuous increase in resolution more shallow convection will also become (partly) resolved and consequently the transport done 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 when this threshold is not met, the model sometimes should resolve at least part of the convection. In these situations the unmodified convection parameterization consumes too much instability, thereby preventing the model from building up convection itself. We therefore need more sophisticated ways to adapt the convection scheme to the gray zone.

[85] presented an approach to make the total turbulent transport, i.e., diffusive plus convective turbulent transport, scale aware. In this approach the parameterized turbulent transport is reduced as a function of a non-dimensionalized grid size. The non-dimensionalized grid size includes the boundary layer height as an estimate of the typical extent of the turbulent transport. [86] showed how this approach could be applied in a convection scheme, and a similar approach is available in HARMONIE-AROME.

Unfortunately, the boundary layer height does not always provide an appropriate indication of the typical scale of the turbulent transport [87]. In practice, this means that the reduction of the parameterized convection according to [85] can be too small, and consequently the model is restricted too much in building up resolved convection itself. Therefore, we are combining the above scale aware convection scheme with another approach, first presented by [88], in which the parameterized convection is shut down if the absolute value of the resolved vertical velocity exceeds a certain value. Here, the threshold on resolved vertical velocity is used as an indication that the model has started to resolve convection by itself. Preliminary results with the combined approach of [85] and [88] in HARMONIE-AROME are promising. The results also suggest that the two updraft framework, as applied in the HARMONIE-AROME convection scheme, provides important refinements in the application of [88].

Note that in the approaches mentioned above we ignored the mass flux assumptions becoming unjustified and the intrinsically stochastical nature of convection, both becoming increasingly important at higher resolutions. To address the stochastic nature of convection in a simple pragmatic way, a stochastically perturbed parameterization ensemble prediction system (SPP-EPS) perturbation can be applied [89].

5.3. Cloud Microphysics

A future addition of the Liquid Ice Multiple Aerosol (LIMA) microphysical scheme [90] is planned. LIMA, the two-moment bulk scheme, is used in other canonical system configurations of the ALADIN-HIRLAM NWP system. It simulates much the same processes and predicts the same hydrometeors as the reference ICE3, but also includes prognostic presentation of the aerosol population. This prognostic aerosol concentration is used in the derivation of the cloud condensation nuclei, cloud droplet number concentration, ice-freezing nuclei, ice crystal, and coated ice-freezing nuclei, which is a cloud condensation nuclei that has frozen. These variables are therefore diagnostic, varying in time. As in the case of ICE3, values for these variables are prescribed. Furthermore, ICE3 does not include coated ice-freezing nuclei in its variables, nor does it include all the microphysical processes which LIMA takes into account, like the Hallett-Mossop process. With that, the addition of LIMA acts as an alternative aerosol option to the near real-time option available in ICE3, and an option for the microphysics in general.

5.4. 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 (so-called ISBA-DIF), and the multi-layer snow scheme (the so-called ES scheme, the Explicit Snow scheme), [40]. We expect a better representation of cooling at night with the multi-layer approach, as was shown in preliminary results by ECMWF and COSMO (personal communication). We also hope for a better representation of the annual cycle of soil moisture, especially during spring. We also plan to apply the Multi-Energy-Balance (MEB) scheme [91,92], 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 from the introduction of the roughness-sublayer (RSL) scheme [93,94], which was recently applied by Shapkalijevski et al. [95]. This scheme considers the canopy-induced drag in surface-atmosphere coupling, as well as the dependency of the displacement height and the roughness length on thermal stability over high vegetation, and thus decreases the positive wind bias over this patch.

Pioneering work has started to produce a new land-cover map with very fine resolution (60 m) over Europe using various sources of information as described in Bessardon et al. [96], Thomas Rieutord [97], Walsh et al. [98], Keany et al. [99]. Bessardon et al. [96] used an agreement-based method to combine information from many maps to overcome variations in semantic, geographical coverage, resolutions, formats, accuracy, and representative periods to produce a new map called ECOCLIMAP-SG+. In addition to ECOCLIMAP-SG+, an agreement score map was generated, which can be interpreted as a quality or uncertainty map. Thomas Rieutord [97] used convolutional neural networks to improve the ECOCLIMAP-SG+ map where the agreement score map showed low quality. In addition, the neural network is able to produce an ensemble of potential land cover maps. Further work is needed to fully test these in HARMONIE-AROME, including technical work on the supplementary physiography datasets of leaf area index, albedo, tree height and lake depths.

5.5. Dynamics and Model Configuration

The horizontal resolution of NWP models has steadily increased over recent decades, with global models now approaching the kilometre scale [100]. Correspondingly, the resolutions of limited-area models will continue to be pushed to sub-kilometre scales, although this effort faces numerous challenges [101].

Model performance at increasing resolution has been an active area of research and development within ACCORD in recent years, with a number of countries now running experimental sub-kilometre 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. The use of a quadratic or cubic spectral grid, discussed earlier, will improve stability by a reduction in the formal spectral resolution, as will off-centring of nonlinear terms in the time discretisation; in particular, the second-order averaging introduced by Simmons and Temperton [102] is available.

More fundamental considerations of the time discretisation scheme may be necessary at extremely high resolutions. As outlined in Section 4, the general ICI approach has been designed to improve stability. More iterations will allow for stable integrations, but at ever increasing costs. So far, testing with the single-step SETTLS-SI approach of HARMONIE-AROME has been successful at grid resolutions of around 500 m. While stable forecasts have been obtained below this, particularly in the case of complex orography and extreme winds this has sometimes required undesirable parameter choices, for example extremely cold values for the ’reference elastic temperature’ in the SI reference state; see Bénard [103] and Bénard et al. [31] for further discussion of this. More investigation is required, and within the wider ACCORD community considerable effort is underway to address these dynamics challenges; for example see Smolíková and Vivoda [104], Burgot et al. [105].

The increase in the horizontal grid resolution desirably has to be followed by a corresponding increase in the vertical grid. Therefore, two options for an increased vertical grid in HARMONIE-AROME have been developed and are being tested. The first option is based on the operational AROME-France 90-level vertical grid (here referred to as $\mathrm{MF}\_90$, Figure 16), 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 (terrain-following grid). The denser representation of the vertical grid (compared to the 65-level grid), especially in the portion of the troposphere closest to the surface, decreases the uncertainty in 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 levels on cloud fraction is shown in Figure 16 top left, where cloud fraction profiles along the blue dashed line shown in the satellite image (top right) are shown for both 65 and 90 (MC_90) vertical levels in HARMONIE-AROME.

However, the increase in vertical resolution (independent from the increase in horizontal resolution) naturally comes with the need for increased computing power for the same domains of interest. This is especially relevant when trying to numerically stabilize and optimize the system at the lowest model level, requiring smaller model time steps. MetCoOp developers have managed to decrease the computing power requirements 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$, Figure 16). This optimization in $\mathrm{MC}\_90$ allows use of the same model time step and iterations in the numerical solver as in the 65-level simulations without negative effects on the numerical stability.4