The aim of this article is to describe the reference configuration of the convection-permitting numerical weather prediction (NWP) model HARMONIE-AROME, which is used for operational short-range weather forecasts in Denmark, Estonia, Finland, Iceland, Ireland, Lithuania, the Netherlands, Norway, Spain, and Sweden. It is developed, maintained, and validated as part of the shared ALADIN–HIRLAM system by a collaboration of 26 countries in Europe and northern Africa on short-range mesoscale NWP. HARMONIE–AROME is based on the model AROME developed within the ALADIN consortium. Along with the joint modeling framework, AROME was implemented and utilized in both northern and southern European conditions by the above listed countries, and this activity has led to extensive updates to the model’s physical parameterizations. In this paper the authors present the differences in model dynamics and physical parameterizations compared with AROME, as well as important configuration choices of the reference, such as lateral boundary conditions, model levels, horizontal resolution, model time step, as well as topography, physiography, and aerosol databases used. Separate documentation will be provided for the atmospheric and surface data-assimilation algorithms and observation types used, as well as a separate description of the ensemble prediction system based on HARMONIE–AROME, which is called HarmonEPS.
There is a strong history of active collaboration between European meteorological institutes on numerical weather prediction (NWP) in order to develop and maintain numerical short-range weather forecasting systems for operational use.
The international research program High Resolution Limited Area Model (HIRLAM) was initiated in 1985 and consists today of the National Meteorological Services (NMSs) from 10 countries: Denmark, Estonia, Finland, Iceland, Ireland, Lithuania, the Netherlands, Norway, Spain, and Sweden, with France as an associate member. Similarly, the collaboration among the NMSs of central Europe, Aire Limitée Adaptation Dynamique Développement International (ALADIN) started in 1991 and consists today of 16 member countries: Algeria, Austria, Belgium, Bulgaria, Croatia, the Czech Republic, France, Hungary, Morocco, Poland, Portugal, Romania, Slovakia, Slovenia, Tunisia, and Turkey.
The ALADIN NWP system is being developed within the frameworks of Action de Recherche Petite Échelle Grande Échelle (ARPEGE) and Integrated Forecasting System (IFS) software, developed jointly by the European Centre for Medium-Range Weather Forecasts (ECMWF) and Météo-France. A more detailed explanation of the ALADIN code architecture and its canonical model configurations, Applications of Research to Operations at Mesoscale (AROME) and Aire Limitee Adaptation/Application de la Recherche a l’Operationnel (ALARO), can be found in P. Termonia et al. (2017, unpublished manuscript).
On 5 December 2005, a cooperation agreement between the ALADIN and HIRLAM consortia was signed with the prime objective “to provide the ALADIN and the HIRLAM Members with a state-of-the-art NWP model for short- and very-short-range forecasting including nowcasting, for both research and development activities and operational usage” (Malcorps and Ågren 2005). In 2014, the ALADIN and HIRLAM consortia further agreed on the formation of a single, united consortium by 2020 and are currently working on this objective. Since 2005, the focus of the HIRLAM research collaboration has been on the convection-permitting scale, and on adapting the AROME model (Seity et al. 2011) for use in the common ALADIN–HIRLAM NWP system, in order to make it accessible for all 26 countries.
The scripting system, which facilitates data assimilation and observation handling, climate generation, lateral boundary coupling, and postprocessing required to run AROME operationally within the HIRLAM countries, is referred to as the HIRLAM–ALADIN Research on Mesoscale Operational NWP in Euromed (HARMONIE) script system. However, the implementation and optimization of AROME for both northern and southern European conditions has led to extensive adaptations and improvements to the model’s physical parameterizations. This was done in order to reduce existing biases and improve the physical description of clouds (mixed phase) and the land surfaces, especially in northern latitude conditions. The model configuration, which uses the updates in the physical parameterizations, has also been referred to as “HARMONIE,” in order to distinguish it from the AROME-France setup. Thus, there is an increased need to clarify and document what is meant by HARMONIE.
The aim of this paper is to describe the reference model configuration of AROME as defined by the HIRLAM consortia: HARMONIE–AROME. It summarizes the changes to the physical parameterizations and dynamics used in HARMONIE–AROME with respect to the description of AROME-France given by Seity et al. (2011) and Brousseau et al. (2016). This paper is limited to the forecast model description of version cycle 40h1.1.1 Separate documentation will be provided for the atmospheric and surface data-assimilation algorithms and observation types used, as well as a separate description of the ensemble prediction system based on HARMONIE–AROME, which is called HarmonEPS.
2. Model dynamics
HARMONIE–AROME uses the same nonhydrostatic (NH) dynamical core as AROME-France, which has been developed by ALADIN (Bubnova et al. 1995; Bénard et al. 2010). It is based on the fully compressible Euler equations (Simmons and Burridge 1981; Laprise 1992). The evolution of the equations is discretized in time and space using a semi-Lagrangian (SL) advection scheme on an A grid and a semi-Implicit (SI) two-time-level scheme, with spectral representation of most prognostic variables based on a double Fourier decomposition. The spectral SI SL scheme originates from the global IFS used operationally at ECMWF (ECMWF 2015a). Horizontal diffusion is applied both by linear spectral diffusion and nonlinear flow dependent diffusion which acts through SL advection and, thus, was given the name semi-Lagrangian horizontal diffusion (SLHD) (Váňa et al. 2008; Bengtsson et al. 2012). Quasi-monotonic operators in the interpolation process are used in order to remove the appearance of negative values for positive definite fields, as well as an unrealistic increase of eddy kinetic energy during the forecast (Seity et al. 2011).
The so-called stable extrapolation two-time-level scheme (SETTLS), specific to the HARMONIE–AROME configuration, is used as the second order two-time-level scheme in order to avoid extrapolation in time of the velocities used for the computation of the trajectories, and for the nonlinear terms, of the evolution equations (Hortal 2002). Furthermore, in order to assure stability of the integrations, a new method for treating the upper-boundary conditions was implemented, using the same Davies–Kallberg relaxation scheme as for the horizontal (Davies 1976). This method makes it possible to use SETTLS also for horizontal resolutions below 1 km, so there is no need to use the so-called predictor–corrector method, which was previously used for very-high-resolution simulations. Another method introduced in order to ensure stability for cases where data assimilation is not used, is to limit the three-dimensional divergence for the first few time steps until equilibrium is reached. This was done in order to avoid occasional crashes when the model starts from model data that is interpolated from a lower-resolution hydrostatic model.
There is also a possibility to run the IFS/ARPEGE software with alternative spectral grids, such as cubic and quadratic grids, and for such grids, a further spectral truncation is performed. Test simulations in HARMONIE–AROME with a cubic grid shows a reduction of the computational time with more than 20%. However, noticeable smoothing and degradation for wind speed in areas of steep orography can be seen compared with the linear grid. Two reasons to consider it nonetheless are 1) to permit running mesoscale ensembles at a reasonable cost, though with somewhat reduced performance with respect to the linear grid, or 2) to reduce model costs of a model upgrade to higher gridpoint resolution. In case of the linear grid, the use of a filter to the vorticity part of the pressure-gradient term is applied in order to eliminate some noise. The method has been applied in the ECMWF model (N. Weidi 2016, personal communication) and is introduced in HARMONIE–AROME with the addition of filtering the pressure departure. In case of quadratic and the cubic grids, the filter is not necessary since the waves modified by the procedure are cut out in those cases.
The Euler equations in AROME-France (and HARMONIE–AROME) are formulated in a terrain-following pressure-based sigma-coordinate system (Simmons and Burridge 1981; Laprise 1992; Bubnova et al. 1995). For the model dynamics, the mean orography may be truncated and smoothed, depending on the transformation between spectral and gridpoint representations. For the linear grid of HARMONIE–AROME, a reduction by a factor of 5 of the shortest wavelength spectrum of the surface elevation is obtained by means of a 16th-order diffusion operator. After smoothing, the atmospheric and surface physical parameterizations refer to the smoothed grid-scale surface elevation. The smoothed or truncated grid-scale surface elevation represents scales somewhat greater than the model’s nominal horizontal resolution.
In the reference cycle of HARMONIE–AROME cycle 40h1.1, lateral boundary conditions are routinely used from the ECMWF model, as opposed to the AROME-France configuration in which the global ARPEGE model is used to provide the lateral boundary conditions (the HARMONIE script system however allows to couple the model to a number of global forecast models). Sixty-five levels are used in the vertical, with model top at ca 10 hPa and lowest level at 12 m. The horizontal resolution is 2.5 km, and the model time step is 75 s.
3. Model physics
The default shortwave (SW) radiation parameterization in AROME-France and HARMONIE–AROME is the Morcrette radiation scheme from ECMWF, IFS cycle 25R1, and contains six spectral intervals (0.185–0.25, 0.25–0.44, 0.44–0.69, 0.69–1.1, 1.1–2.38, and 2.38–4.00 μm). The default longwave (LW) radiation scheme contains 16 spectral bands between 3.33 and 1000 μm. This uses the Rapid Radiative Transfer Model (RRTM) of Mlawer et al. (1997). Both the SW and LW schemes are described in the IFS (ECMWF 2015b) and the mesoscale research model Meso-NH (Mascart and Bougeault 2011) documentation. Because of computational constraints the full radiation calculations are currently performed every 15 min. The more affordable single-band radiation schemes from ALARO physics (ACRANEB2; Mas̆ek et al. 2016; Geleyn et al. 2017) and HIRLAM (HLRADIA; Savijärvi 1990; Wyser et al. 1999), which can be run at each time step, have also been implemented in HARMONIE–AROME for experimentation purposes. Hereafter, unless stated otherwise, when we refer to the SW and LW radiation schemes the default parameterizations are implied. The following description focuses on SW parameterizations because the LW RRTM scheme is applied with minimal modifications.
The clear-sky SW radiative transfer is calculated using the Fouquart and Bonnel (1980) two-stream equations. The reflectance, absorption, and transmittance of the clear-sky fraction of the atmospheric layers are calculated in a similar manner to that outlined in Coakley and Chylek (1975). The cloudy-sky SW computations are done using the delta-Eddington approximation of Joseph et al. (1976). The radiative transfer calculations use the inherent optical properties [(IOPs): optical thickness, single scattering albedo (SSA), asymmetry factor (g)] of cloud particles (prognostic specific cloud liquid and cloud ice content), aerosols (monthly climatologies), and atmospheric gases (prognostic H2O, a fixed composition mixture of CO2, N2O, CH4, and O2, monthly climatologies of O3).
A variety of options for the parameterization of cloud particle size and shape, and for the consequent derivation of cloud optical properties, are available within the IFS radiation scheme. Table 1 shows the choices recommended for the HARMONIE–AROME reference cycle 40h1.1. These choices differ from the defaults used in the AROME-France cycle 40t1 setup (Seity et al. 2011). In particular, we have introduced an improved cloud liquid optical property scheme (Nielsen et al. 2014; Gleeson et al. 2015), which is based on detailed Mie theory computations. In comparison with the accurate one-dimensional radiative transfer model, Discrete Ordinate Radiative Transfer model (DISORT) (Stamnes et al. 1988, 2000), the new cloud liquid optical property scheme is shown superior to the previous scheme. Furthermore, in the previous version of the model, the cloud inhomogeneity factor was assumed to be 0.7 in order to account for a variability of cloud in a grid box. This assumption is no longer valid with increased model grid resolution (R. Hogan 2014, personal communication; A. M. Townsend 2015, unpublished manuscript); thus, we have assumed that the clouds are homogeneous where present in a grid cell. As a first approximation, the radiative effect of precipitating graupel and snow particles is included by assuming these to have the same inherent optical properties as cloud ice; this was done in conjunction with the inclusion of the cloud microphysics updates described in section 3b. Both the Nielsen scheme and the reduced cloud inhomogeneity lead to a decrease in the downwelling SW radiation flux (Gleeson et al. 2015). However, an overestimation of low-level clouds have been reduced in the new cycle by introducing stronger mixing in the boundary layer, using a new turbulence scheme (described further down); thus, the overall impact of the radiation updates in the new cycle are rather neutral, albeit more correct. In each column, a maximum-random cloud overlap is assumed in the vertical.
The direct SW radiative effect of aerosols is calculated using vertically integrated aerosol optical depth (AOD) at a wavelength of 550 nm (AOD550) and the following aerosol IOPs: AOD spectral scaling coefficients and spectral single scattering albedo (SSA) and asymmetry factor (g). The indirect radiative effect of aerosols due to cloud particle formation is not included in the current version of HARMONIE–AROME. Monthly climatologies of AOD550 of land, sea, desert, and urban tropospheric aerosols from the Tegen et al. (1997) climatology are used along with background stratospheric aerosols in a similar manner to the IFS model (ECMWF 2015b). These are distributed among the model levels using the Tanré et al. (1984) climatological vertical profiles for each aerosol type (see Gleeson et al. 2016; Toll et al. 2016). The spectral dependence of AOD, SSA, and g for each aerosol type is parameterized following Hess et al. (1998). Toll et al. (2016) showed that in the Tegen et al. (1997) climatology the AODs are underestimated over Europe compared to more recent datasets, especially near the Atlantic Ocean coasts. This leads to some overestimation of the clear-sky SW irradiance at the surface.
Grid-scale surface SW albedo and LW emissivity, required as a boundary condition by the radiation parameterizations, are based on surface characteristics given by the 1-km-resolution ECOCLIMAP database (Faroux et al. 2013) and processed by the Surface Externalisée (SURFEX) externalized surface scheme (Masson et al. 2013; described in more detail in section 3e). Ultraviolet, visible, and near-infrared values of the surface albedo are mapped to the six SW spectral bands for the ECMWF IFS scheme (IFS cycle 25R1; ECMWF 2015b) though the UV albedo is unused in practice. The single-band ACRANEB2 scheme is interfaced to SURFEX using one SW spectral interval on both sides. By default, SURFEX assumes the same value for the direct and diffuse albedo for each band. We improved this by applying an empirical correction, which depends on the solar zenith angle (SZA), to the diffuse albedo (αdif) in order to derive the direct beam albedo (αdir): αdir = αdif + 0.2/[1 + cos(SZA)] − 0.12. This correction was imported from the HIRLAM model (Unden 2002).
Diagnostic output from the radiation parameterizations includes accumulated spectrally averaged downwelling SW global, direct and direct normal irradiances at the surface, and net SW and LW radiative fluxes at the top of the atmosphere, at the surface, and on each model level. The downwelling diffuse SW radiation can be obtained from the difference between the global and direct radiation at the surface. Diffuse radiation includes both cloudy- and clear-sky contributions, whereas a small part of the direct radiation is assumed to come from the cloudy sky. In the original IFS cycle 25R1 scheme, direct and clear-sky radiation were assumed to be identical as were diffuse and cloudy-sky radiation fluxes, which are incorrect assumptions.
b. Clouds and cloud microphysics
The microphysics scheme used in AROME-France and HARMONIE–AROME is a one-moment bulk scheme, which uses a three-class ice parameterization, referred to as ICE3, originally developed for Meso-NH (Pinty and Jabouille 1998; Lascaux et al. 2006). It contains the following solid hydrometeors as prognostic variables: cloud ice, snow, and a combination of graupel and hail. Graupel and hail may be separated by using an own prognostic variable for hail, but this is still in research mode at present. The other prognostic variables used in the cloud microphysics scheme are water vapor, cloud liquid water, and rain. All hydrometeors are advected horizontally by a semi-Lagrangian scheme and vertically by a sedimentation scheme described in detail in Bouteloup et al. (2011). Three-dimensional cloud fraction is not a prognostic variable but instead is determined using a statistical cloud and condensation scheme (Bougeault 1982; Bechtold et al. 1995).
Some weaknesses in the original scheme have been detected, particularly in the stable boundary layer during winter over northern Europe. In these situations, the model generates ice too quickly when temperatures in the clouds are between −5° and −10°C, where some supercooled liquid would be expected. Furthermore, at temperatures lower than −20°C, spurious clouds are often present at the lowest model level and may be treated as “fog” by users of the model output, while observations show clear skies. The reason for this is that the original scheme, by construction, removes most supersaturation with respect to ice in regions where the temperature is below −20°C and forms ice clouds, while in reality, supersaturation in such conditions is common, since ice clouds are formed at a much slower rate than the typical time step used in the model.
To address these weaknesses, substantial updates have been made to the cloud microphysics scheme under the option “OCND2,” which was introduced in order to improve clouds in cold conditions, described in more detail in Müller et al. (2017). The main difference compared to the original scheme is that in OCND2, the fraction of the grid box with cloud ice (and with supersaturation with respect to ice) is no longer handled by the large-scale condensation and thermodynamic adjustment scheme; instead, it is parameterized using a cloud scheme based on the critical relative humidity with respect to ice. In the original scheme, even though cloud ice is a prognostic variable, it is treated similarly to a diagnostic quantity as it is a function of temperature only. In OCND2, only cloud water is handled by the large-scale condensation and thermodynamic adjustment scheme; cloud ice is treated by the rest of the ICE3 microphysics, which includes sublimation, evaporation, and interactions with other water species.
Besides this improved separation between the fast liquid processes and the slower ice water processes, some other updates are included in the OCND2 scheme: a reduction of the deposition rate of the ice-phase water species, a correction of the total cloud cover to address the lower optical thickness of ice clouds compared with water clouds, and a reduction of the ice nucleolus concentration in temperatures between 0° and −25°C. Furthermore, the process of rain drop activation from cloud droplets (autoconversion) is parameterized using the “Kogan autoconversion” parameterization (Khairoutdinov and Kogan 2000) as opposed to the Kessler (1969) scheme used in AROME-France.
An example of the impact of OCND2 on modeled cloud liquid and ice phases can be seen in Fig. 1. The cloud liquid and cloud ice phases from model runs with and without OCND2 are compared with observations from the Hyytiälä station located close to Helsinki, Finland, for the month of February 2014. Hyytiälä is an “ARM mobile facility,” and data from the site are used within the CloudNET project (Illingworth et al. 2007) to evaluate the representation of clouds in climate and weather forecast models. The observed cloud liquid water content is calculated within CloudNET using both cloud radar and lidar as well as dual-wavelength microwave radiometers. More information about the CloudNET method can be found in Illingworth et al. (2007) and references therein. The OCND2 scheme improves the representation of mixed-phase and pure ice clouds in the model in the wintertime by increasing the amount of liquid water in low-level clouds in cold conditions and decreasing the amount of ice water content (ice + graupel + snow) in low-level clouds, such that they are closer to the observed (by CloudNET) values of liquid and ice water content (Fig. 1). In the summertime, the impact is not as large; however, more supercooled liquid can be seen in higher-altitude clouds with the OCND2 scheme in the summertime, which is closer to the observations (not shown).
The new representation of mixed-phase clouds has led to an improvement in many of the meteorological fields. In particular it has led to a reduction of a cold bias in the 2-m temperature during wintertime (over Scandinavia) and reduction of an existing dry bias in relative humidity throughout the lower atmosphere in winter (see Fig. 2).
In earlier versions of HARMONIE–AROME a number of persistent deficiencies in the representation of the boundary layer could be observed: too low boundary layer heights and clouds base, too much cloud cover, and too much fog, in particular over sea (de Rooy 2014). One-dimensional versions of AROME-France and HARMONIE–AROME participated in the Atlantic Stratocumulus to Cumulus Transition Experiment (ASTEX) intercomparison study. Results of HARMONIE–AROME with the CBR turbulence scheme used in AROME-France (Cuxart et al. 2000; Seity et al. 2011) revealed a substantial underestimation of the cloud-top entrainment by the turbulence scheme for this case (de Rooy 2014).
Therefore, a new turbulence scheme, HARMONIE with RACMO Turbulence (HARATU), which has a larger cloud-top entrainment, has been implemented in HARMONIE–AROME cycle 40h1.1. HARATU is based on a scheme that was originally developed for use in the regional climate model RACMO (van Meijgaard et al. 2012; Lenderink and Holtslag 2004). Similar to the CBR scheme, it uses a framework with a prognostic equation for the turbulent kinetic energy (TKE) combined with a diagnostic length scale. The TKE equation includes source (+) and sink (−) terms due to wind shear (+), buoyancy (+ for unstable and − for stable conditions), transport (locally + or −, but no net effect), and dissipation of TKE (−).
Compared to the CBR scheme, there are considerable changes in the length-scale formulation and the constants used. In the CBR scheme there is one “master” length scale, which is multiplied by a number of stability dependent functions. In HARATU the stability corrections are part of the length-scale formulation. As such, there are different length scales for heat and momentum. Also, the numerical implementation of the TKE equation has been changed from “full” levels (where the temperature, moisture and wind are computed) to “half” levels (where the fluxes are computed). This choice avoids unnecessary vertical interpolations in the computation of the turbulent fluxes and the source and sink terms of the TKE equation. In particular, in the case of strong gradient, this gives more reliable estimates of the turbulent fluxes and vastly improves the cloud-top entrainment.
The length-scale formulation in HARATU essentially consists of two parts: one for stable condition and one for near-neutral to convective conditions [see Lenderink and Holtslag (2004) for an extensive description]. The stable length-scale formulation is the commonly used buoyancy-based length scale given by the square root of TKE divided by the vertical stability (Deardorff 1980; Baas et al. 2008). The neutral–unstable length scale consists of vertical integrals of stability dependent functions. This is done in an upward and a downward computation, and the resulting upward and downward length scale are averaged to obtain the neutral/unstable length scale. The stability functions use the Richardson number (Ri), which allows us to match with surface layer similarity for near-neutral conditions (Lenderink and Holtslag 2004).
Stability coefficients take into account moist processes; that is, the effect of latent heat on stability due to condensation or evaporation of cloud droplets. In general, these moist processes introduce a strong coupling between the turbulence scheme and the cloud and condensation scheme, and this makes moist turbulence schemes very susceptible to numerical instability and noise (Lenderink et al. 2004; Lenderink and van Meijgaard 2001). The length-scale formulation here is rather insensitive to those numerical instabilities mainly because of its formulation where the integral over stability is used. Having this formulation generally produces smooth and continuous results, in particular in the presence of clouds (Lenderink and Holtslag 2004; Lenderink et al. 2004).
With respect to the original turbulence scheme described in Lenderink and Holtslag (2004), a few important modifications have been made when implementing the scheme in HARMONIE–AROME. To combine the scheme with the dual mass-flux scheme described below, the stability functions for the near-neutral/unstable length-scale formulation had to be modified for the following reason. The mixing of heat due to the mass-flux scheme in a convective boundary layer leads to a slightly stable temperature profile in the upper part of the mixed layer, which is consistent with large-eddy simulation (LES) of a convective boundary layer (de Roode et al. 2004). With the original formulation, this leads to a strong and unrealistic reduction of the mixing length in the upper part of the convective boundary layer. To avoid this, we adjusted the stability functions using a first-order approximation of the change in the profile due to the mass-flux contribution. Also, a small modification was made to avoid the discontinuity in the Richardson number in the case of vanishing wind shear.
The maximum wind speeds over land in strong wind conditions (>10 m s−1) with the original formulation described in Lenderink and Holtslag (2004) turned out to be approximately 10% lower than with the CBR scheme and appeared to be too low compared to measurements. For this reason, we performed a small retuning of the scheme by enhancing the mixing length near the surface (effectively by 20% for neutral conditions) and adjusting the “downward” length-scale formulation, leading to a more effective downward mixing of momentum in the case of strong winds over land. This modification has almost no influence over sea, and for weak-to-moderate (below 10 m s−1) wind speeds over land.
The implementation of HARATU has considerably reduced the cloud cover and resulted in an increase in clouds’ base height compared to the CBR scheme (de Rooy and de Vries 2017). Furthermore, the HARATU scheme also considerably improves the wind climatology of the model. As an example, over the Netherlands domain, the bias in the diurnal cycle of the mean wind speed is almost zero during the whole day, whereas the previous turbulence scheme as well as the scheme used in the operational ECMWF model (cy41, T1279) shows a clear diurnal signal in the bias over the same domain (see Fig. 3a). This reduction of 10-m wind speed bias has also been seen in all of the other domains running HARMONIE–AROME operationally (not shown). The scheme also improves the standard error and absolute error of the wind forecast (Fig. 3b). In addition, the wind shear in the lower boundary layer is better captured compared to tower observations from the Cabauw measurement site (see Figs. 3c,d). Further evaluation of the wind speed over sea using scatterometer data shows clear improvement to the CBR scheme (de Rooy and de Vries 2017).
At 2.5-km resolution, deep convection is expected to be roughly resolved and explicitly represented by the model’s nonhydrostatic dynamics; thus, in HARMONIE–AROME there is no parameterization of deep convection. However, shallow convection still needs to be parameterized. For this, usually a mass-flux framework is applied consisting of one or more updrafts, which transport heat, moisture, and momentum. The convective transport is proportional to the difference between the updraft properties and the environment, times the amount of mass transported by the updraft. The upward mass flux M is described by a simple budget equation:
where ε is the fractional entrainment, describing the inflow of environmental air into the updraft (herewith diluting the updraft) and δ is the fractional detrainment describing the outflow of updraft air into the environment. These coefficients can be considered as the key elements in a convection scheme. HARMONIE–AROME uses a different scheme for shallow convection than AROME-France (Seity et al. 2011), called EDMFm. According to the original ideas of Siebesma and Teixeira (2000), Soares et al. (2004), Siebesma et al. (2007), and Rio and Hourdin (2008), the mass-flux concept can be applied in a so-called eddy diffusivity mass-flux (EDMF) framework. The eddy diffusivity, or turbulence part is given by the turbulence scheme HARATU explained above. Here we discuss the mass-flux component, which describes the transport by cloudy as well as dry (unsaturated) updrafts. The focus will be on what distinguishes EDMFm in HARMONIE–AROME from the mass-flux scheme used in AROME-France, which is referred to as EDKF (Pergaud et al. 2009).
1) Dual mass flux
Contrary to EDKF, EDMFm uses a dual mass-flux approach in which two updrafts are distinguished: a dry updraft that never reaches the lifting condensation level and a moist updraft that condenses and becomes a cloud. As schematically illustrated in Fig. 4, three different convective boundary layer regimes are considered. In contrast to Neggers et al. (2009) where the subdivision between dry and moist updraft fractions is flexible, EDMFm uses fixed values only depending on the regime.
The scheme starts with the initialization of the excess of the updrafts (Neggers et al. 2009). Subsequently, a vertical velocity equation is used to determine updraft vertical velocity and the corresponding height to which the updraft can penetrate. This provides the inversion, or cloud-base height zi, and cloud-top height zt. The applied formulation of the vertical velocity equation [based on Siebesma et al. (2007), Simpson and Wiggert (1969), and de Rooy and Siebesma (2010)] is recently also supported by de Roode et al. (2012).
With the inversion height known, the profiles of entrainment rate are defined because they are functions of just z and zi (see Fig. 4). For the dry updraft we adopt the ε formulation of Siebesma et al. (2007), based on LES results for the dry convective boundary layer. The ε profile for the moist updraft in the subcloud layer extrapolates the work of Siebesma et al. (2007) for dry updrafts that stop at zi, to the larger subcloud thermals that do not stop at zi. These larger, faster rising thermals are associated with smaller entrainment rates in comparison with the dry updraft (see Fig. 4b). For the moist updraft, the value of entrainment at cloud base scales with as suggested by LES, reflecting that deeper mixed layers can be associated with higher vertical velocities and larger thermals (de Rooy and Siebesma 2010).
Finally, the entrainment formulation in the cloud layer decreases with height as z−1 (Siebesma et al. 2003; de Rooy and Siebesma 2008). A small refinement in EDMFm concerns the value of ε at cloud base, which is connected to the value of ε at the top of the subcloud layer, and thus depends on the mixed layer height.
2) Parameterization of detrainment
The essential difference between EDMFm and EDKF concerns the parameterization of the detrainment in the cloud layer. As first pointed out by de Rooy and Siebesma (2008), variations in the mass-flux profile from case to case and hour to hour can be almost exclusively related to the fractional detrainment (δ). This is supported by numerous LES studies, revealing orders of magnitude larger variations in δ than in ε (e.g., Jonker et al. 2006; Derbyshire et al. 2011; Böing et al. 2012; de Rooy et al. 2013). Apart from this empirical evidence, the much larger variation in δ and its strong link to the mass flux is explained by theoretical considerations in de Rooy and Siebesma (2010). For the first time, the implications of the aforementioned considerations are used in an operational scheme. As a result, EDMFm behaves fundamentally different than other operational schemes. For example, a Kain–Fritsch type scheme (Kain and Fritsch 1990; Kain 2004), like the EDKF option used in AROME (Pergaud et al. 2009), in which ε and δ vary in an opposite but similar manner to environmental conditions, is not able to capture the order of magnitude variations in δ due to cloud layer depth (de Rooy et al. 2013). In EDMFm this cloud layer depth dependence is included by considering the mass-flux profile in a nondimensionalized way (de Rooy and Siebesma 2008). Apart from the cloud layer depth, the detrainment value is influenced by environmental conditions. For this dependency, we use a parameter called χc (Kain and Fritsch 1990). As shown by de Rooy and Siebesma (2008), χc increases with the relative humidity and buoyancy excess of the updraft. Therefore, the LES-based functional dependence of δ (or mass-flux profile) on χc is physically plausible: high values of χc can be associated with large clouds, with high updraft velocities that have large buoyancy excess and/or clouds rising in a humid environment. Accordingly, high χc values correspond to small detrainment values and slowly decreasing mass flux with height. Further details on the detrainment formulation can be found in de Rooy and Siebesma (2008).
Describing lateral mixing in the cloud layer in this way is supported by observational [see, e.g., Lamer et al. (2015)] and LES studies. The most convincing support from LES can be found in Böing et al. (2012), who used 90 LES runs to explore the sensitivity of ε, δ, and the mass-flux profile in deep convection to a broad spectrum of relative humidities and stability of the environment. This study confirms the much larger variation of δ and its strong link with the mass-flux profile. Figure 5, from Böing et al. (2012), shows the results of parameterizing the mass-flux profile according to de Rooy and Siebesma (2008). Another difference between HARMONIE–AROME cycle 40h1.1 and AROME-France is the method by which convection influences the total cloud cover and subgrid liquid and ice water content. In HARMONIE–AROME the method proposed by Soares et al. (2004) is applied, in which the mixing from turbulence and convection is used to produce the variance of the distance to saturation in the statistical cloud scheme. In AROME-France only turbulence contributes in this way to the variance, whereas the impact of convection on the cloud cover and subgrid liquid and ice water content is assumed to be proportional to the updraft area fraction of the mass flux of the previous time step.
The surface physics in AROME-France and HARMONIE–AROME is simulated by the surface scheme named SURFEX (Masson et al. 2013), which is a surface modeling platform developed mainly by Météo-France in cooperation with the scientific community. SURFEX is composed of various physical models for natural land surface, urbanized areas, lakes, and oceans. It also simulates chemistry and aerosols surface processes and can be used for assimilation of surface and near-surface variables. SURFEX has its own initialization procedures and can be used in standalone mode and coupled to an atmospheric model (Masson et al. 2013). In SURFEX, each model grid box is represented by four surface tiles: sea or ocean, lakes, urban areas, and nature (soil and vegetation). The nature tile can further be divided into several so-called patches depending on vegetation type. Each surface tile is modeled with a specific surface model and the total flux of the grid box results from the addition of the individual fluxes weighted by their respective fraction.
HARMONIE–AROME cycle 40h1.1 uses SURFEX version 7.3. The exchange of energy and water between the land surface and the atmosphere above is simulated by the Interactions between Soil, Biosphere, and Atmosphere (ISBA) scheme with a force–restore approach (Boone et al. 1999) in combination with the Douville snow scheme Douville et al. (1995). Currently, the surface characteristics are aggregated, and only one patch is used in the flux calculations on the nature tile. For the sea tile, the Exchange Coefficients from Unified Multi-campaign Estimates (ECUME) scheme by Belamari (2005) is used over water and the “Simple Ice Model” for sea ice, as described below, has been added in the latest version of HARMONIE–AROME. For the inland water tile (lakes and rivers) the Charnock (1955) formula is used over water. The lake surface temperature is initialized by deep soil temperature (extrapolated if necessary) and is kept constant during the forecast. For water temperature below the freezing point, surface properties for snow is applied; for example, surface momentum roughness is set to 0.001 m and albedo to 0.85. Finally, the urban tile is simulated by the Town Energy Balance (TEB) model (Masson 2000).
Over the water and sea tiles, diagnostic quantities at 2 and 10 m are calculated by interpolating between atmospheric forcing variables and surface temperature and humidity variables. Over land, the surface boundary layer (SBL) scheme by Masson and Seity (2009) is used. The 1D prognostic turbulence scheme calculates TKE, wind, temperature, and humidity on six vertical levels 0.5, 2, 4, 6.5, 9, and 12 m above ground. The motivation for using the SBL scheme is to improve the performance in stable situations. However, experiences show that while the scheme might give realistic and low temperatures in some situations, it can also yield much-too-low temperatures in some situations, often in combination with too-weak winds. The physiography databases related to land use, topography, and clay/sand are currently all revisited for domains used within the HIRLAM consortia. The default surface land-cover physiography in cycle 40h1.1 is based on ECOCLIMAPv2.2 (Faroux et al. 2013). Modifications have been included in places where the ECOCLIMAP description was found to be suboptimal (e.g., over the permanent snow areas over Norway). Also, over Greenland and Iceland, local modifications of these databases (e.g., permanent snow areas and leaf-area index) have shown significantly improved results in near-surface wind and temperature scores, demonstrating the importance of carefully looking into the physiography used. The default surface topography is based on Global Multi-resolution Terrain Elevation Data (GMTED2010) (Danielson and Gesch 2011). The default clay and sand proportions are still based on FAO (FAO 2006) since the newer Harmonized World Soil Database (HWSD) (Nachtergaele et al. 2012) shows dubious values over Scandinavia. A monthly climatologies of vertically integrated optical depth of four aerosol species (Tegen et al. 1997) dataset is introduced to the forecast model along with the physiography and topography data.
Previous versions of HARMONIE–AROME treated areas covered by sea ice in a quite simplified manner, using a constant value for sea ice surface temperature during the whole forecast. However, it was found that such a configuration led to a noticeable bias of 2-m temperature over ice covered areas which grew with increasing forecast lead time. To solve this problem, the Simple ICE model (SICE) was introduced in HARMONIE–AROME cycle 40h1.1 (not activated by default, but switched on and used by some services where sea ice is a big part of the domain). SICE is built on top of SURFEX’s soil heat-diffusion solver and represents a layer of sea ice with fixed thickness and prognostic temperature within the ice slab. The ice pack is divided into a number of layers in order to solve the heat diffusion. Here the uppermost layer is defined by the ice surface temperature, which is derived from the thermal balance equation, and the lowermost layer holds the freezing point temperature. Ice covered areas are determined by the ice concentration field provided by an external source. The flux from a sea tile grid cell is calculated as the weighted contribution from the ice and open water schemes. HARMONIE–AROME cycle 40h1.1 uses SICE configuration with 0.75-m-thick ice slab divided into four layers.
Figure 6 shows the impact of the SICE scheme over seven stations in the Gulf of Bothnia for the time period 1–31 March 2013, forecasts initialized at 0000 UTC. Black curves represent the HARMONIE–AROME model without the SICE scheme, and red curves are with the SICE scheme. Here it can be seen that the SICE scheme improves the forecast of surface pressure, temperature, and wind speed. Temperature is particularly improved, with a better daily cycle and a smaller model bias.
4. Future developments in HARMONIE–AROME
In the dynamics, a development for the near future is the introduction of vertical finite elements, which has been done in close collaboration with the ALADIN consortia. An advantage of this vertical discretization is the exclusive use of full levels, skipping a computational mode created from the interpolation from full levels to half levels. The vertical finite elements were successfully introduced by Untch and Hortal (2004) in the hydrostatic ECMWF IFS model. The implementation of vertical finite elements has also been extended to work for the nonhydrostatic dynamics. However, the nonhydrostatic model needs to solve a constraint called C1 involving the SI vertical operators, which is not present in the hydrostatic version. If this constraint is fulfilled then one can write the SI set of linear equations only in terms of the vertical divergence. This has been done in the finite difference discretization, but unfortunately, the C1 constraint is not guaranteed in the construction of the finite element operators. An iterative method was developed by (Vivoda and Smolíková 2013) in order to relax the C1 constraint. Based on this approach a theoretical development to solve the C1 constraint was proposed by Subias (2015), which has been tested in HARMONIE–AROME in cycle 40h1.1. The implementation of the vertical integral operators involved in the C1 constraint shows no impact in 3D tests when they replace the default operators in a purely finite difference configuration. Thus, these operators are good candidates for use in a nonhydrostatic finite element configuration because they satisfy the C1 constraint. However, the tests are very sensitive to the choice of the vertical levels, so future work is needed to adapt the scheme to work for any given set of levels.
We will also continue to seek adaptations of the semi-Lagrangian method that conserve mass better but do not involve a large increase in computational cost. In this context we will explore further the use of the Continuous Mapping about Departure points (COMAD) scheme (Malardel and Ricard 2015), which introduces a correction applied to the standard interpolation weights in the SL scheme and takes into account the deformation of the air parcels along each direction of interpolation. The scheme is already used in AROME-France, and recent tests in HARMONIE–AROME show that the scheme helps to reduce excessive buildup of cloud hydrometeors in isolated grid points. As a next step, it may be considered to enhance mass conservation for individual atmospheric components in the SL treatment of the mixing ratio equations.
On a longer time scale, the semi-implicit time stepping scheme will be reconsidered at very high (subkilometer) resolutions; here we will look at steep slope behavior and computational performance and assess the potential of the alternative horizontal explicit vertical implicit (HEVI) scheme (Lock et al. 2014).
The ACRANEB2 radiation scheme from the ALARO model (Mas̆ek et al. 2016; Geleyn et al. 2017) and the HLRADIA radiation scheme from the HIRLAM model (Savijärvi 1990; see also Nielsen et al. 2014) are now available for testing in HARMONIE–AROME; we will investigate whether the ability to have fast interactions between clouds and radiation and the surface and radiation are of greater importance for model performance than accounting for the spectral details of clear-sky radiation. Also, it will be investigated whether these radiation schemes can be used in individual ensemble members for HARMONIE–AROME ensemble simulations.
For direct and indirect aerosol parameterizations, input information on the atmospheric aerosol distribution and its optical and chemical properties is required. To improve the present aerosol climatology (Tegen et al. 1997), an update of the aerosol climatology to the Monitoring Atmospheric Composition and Climate (MACC) reanalysis (Inness et al. 2013) dataset, which includes assimilated AOD measurements, is considered, as well as the use of real-time aerosol data from Copernicus Atmosphere Monitoring Service (CAMS). Updated aerosol data will be utilized both by the radiation and cloud-precipitation microphysics parameterizations.
c. Clouds and microphysics
Improvements in the parameterizations of cloud microphysics and hydrometeor interactions within the clouds are sought with the aim to enhance forecast accuracy for extreme precipitation events and the prediction of fog and low clouds. In future model versions, we will explore the two-moment cloud microphysics scheme Liquid Ice Multiple Aerosols (LIMA) developed by the Meso-NH community and at Météo-France (Vié et al. 2016), which has been derived from the ICE3 microphysics scheme. The scheme introduces prognostic variables for droplet number concentration for cloud, rain, and ice and allows for a more realistic description of cloud–aerosol interactions. We will work in close collaboration with the ALADIN consortia on testing the scheme with various sources of aerosol, MACC reanalysis, and real-time aerosol analysis from Copernicus.
d. Turbulence and convection
Work to prepare the model for operational use at increased resolution (100 layers, 0.5–1.3 km) will be a priority the coming years, as well as exploring the model behavior in the gray zone of shallow convection and turbulence. Experiences of using the model at these resolutions can be drawn from, for example, Brousseau et al. (2016) and Honnert et al. (2011). There are also plans to understand better the behavior of the HARATU scheme in the stable boundary layer and whether the model can be improved further in this regime.
A number of deficiencies in the performance of the HARMONIE–AROME configuration can be attributed to surface processes and physiography issues. These concern for instance a cold and humid spring bias over northern Europe, a humid and cold early spring in western Europe, followed by a warm and dry late spring and summer period and a shift in temperature climate nearby deep and large lakes.
One promising step toward an improved surface description is to increase the number of patches over the nature part from one to two; that is, subdivide the nature tile into a forest and an open land patch, respectively. With one patch all surface properties are averaged to land-averaged values while with two patches each patch is given its unique surface properties. Preliminary results look promising and show an increase in Bowen ratio and a reduced problem with respect to the too-humid spring conditions over northern Europe but also show a reduced winter temperature bias over southwestern Europe.
However, achievement of a more complete solution requires utilization of more advanced surface modules that have become available in SURFEXv8 in cycle 43 of the ALADIN–HIRLAM NWP system. These processes concern multilayer soil and snow schemes and an explicit canopy treatment where the canopy vegetation, energy-budget-wise, is separated from the soil and snow beneath. Assessing the potential of these schemes should be done in close connection to the corresponding data-assimilation methods. For their initialization, the surface data assimilation needs to be based on more sophisticated algorithms than the present Optimum Interpolation, such as the extended Kalman filter (EKF). This will also make it possible to utilize a wide range of remote sensing products for surface data assimilation.
The temperature problem connected to lakes will be addressed by activating the lake model FLake (Mironov et al. 2010) for all inland water (lakes and rivers). FLake is already used operationally in NWP by, for example, COSMO (Mironov et al. 2010) and ECMWF (ECMWF 2015b), and in climate applications of the ALADIN–HIRLAM system using the HARMONIE Climate configuration (Lind et al. 2016). The most important physiography information for lakes is the lake depth. For this we use the Global Lake Depth Database (GLDB) by Choulga et al. (2014). The large thermal inertia of lakes does also require a careful initialization. For this, the global lake climatology by Kourzeneva et al. (2012) is used. In the coming years, an EKF lake data-assimilation scheme will be developed in which satellite and in situ observations of lake surface temperature and ice cover can be assimilated.
Furthermore, parameterizations of the orography impact on surface-level radiation fluxes based on Senkova et al. (2007) [see also Rontu et al. (2016)] have been prepared within SURFEX. Coupling of these parameterizations to the full HARMONIE–AROME model will be tested within the next cycle of the system.
f. Coupling with sea surface and ocean
The use of the HARMONIE–AROME model is gradually extending beyond NWP to include more Earth system components and to longer time frames. Preliminary experiments in Norway have indicated that coupling HARMONIE–AROME with the wave model WAM is beneficial in the sense that it reduces the systematic increase in near-surface wind bias for strong winds, which has been observed in verification against scatterometer data and buoys (Süld et al. 2015). For this reason, and considering the relevance of this benefit for accurately predicting polar lows, a two-way coupling with WAM has been incorporated in model simulations over an Arctic domain used in operations by Met-Norway. Having a two-way coupled ocean–atmosphere model is a possible step for the future, in particularly for the regional climate modeling community. In this regard, ALADIN partners, in particular Météo-France and Croatia, have made progress: at Météo-France by coupling AROME-France to the ocean model NEMO (Madec et al. 2015), using the OASIS coupler (Valcke 2013), which exists within SURFEX (Masson et al. 2013)—this work is described in the Ph.D. thesis of Rainaud (2015); in Croatia by testing a two-way coupling between the ALADIN model on the atmosphere side, and an Adriatic setup of Princeton ocean mode, POM (Blumberg and Mellor 1987), on the ocean side as described by Ličer et al. (2016). Based on these experiences we aim to make progress on atmosphere–ocean coupling in the coming years in the ALADIN–HIRLAM NWP system.
We thank our colleagues within the HIRLAM and ALADIN consortia for active and stimulating collaboration leading up to the completion of this manuscript.
The ALADIN–HIRLAM limited area system is part of the code base of the IFS/ARPEGE system of ECMWF and Météo-France. The cycle number refers to the main model version as released by ECMWF. The subsequent letter (r, t, or h) refers to new model updates released by ECMWF, ALADIN, or HIRLAM consortia, respectively.