Enviro-HIRLAM online integrated meteorology-chemistry modelling system strategy, methodology, developments and applications (v7.2)

. The Environment – High Resolution Limited Area Model (Enviro-HIRLAM) is developed as a fully online integrated numerical weather prediction (NWP) and atmospheric chemical transport (ACT) model for research and forecasting of joint meteorological, chemical and biological weather. The integrated modelling system is developed by the Danish Meteorological Institute (DMI) in collaboration with several European universities. It is the baseline system in the HIRLAM Chemical Branch and used in several countries and different applications. The development was initiated at DMI more than 15 years ago. The model is based on the HIRLAM NWP model with online integrated pollutant


Introduction
During the last decades, a new field of atmospheric modelling -the chemical weather forecasting (CWF) -has been quickly developing and growing.However, in most of the current studies, this field is still considered in a simplified concept of the offline running of atmospheric chemical trans-port (ACT) models with operational numerical weather prediction (NWP) data as a driver (Lawrence at al., 2005).A new concept and methodology considering the "chemical weather" as two-way interacting nonlinear meteorological and chemical/aerosol dynamics processes of the atmosphere have been recently suggested (Grell et al., 2005;Baklanov and Korsholm, 2008;Baklanov, 2010;Grell and Baklanov, 2011).First attempts at building online coupled meteorology and air pollution models for environmental applications were done in the 1980s; see Baklanov (1988), Schlünzen and Pahl (1992), Jacobson (1994).For climate applications, the first coupled chemistry-climate models were developed and used in the 1990s; cf.Jacobson (1999Jacobson ( , 2002)), de Grandpré et al. (2000), Steil et al. (2003), Austin and Butchart (2003).A more detailed overview of the history and current experience in the online integrated meteorology-chemistry modelling, the importance of different chains of feedback mechanisms for meteorological and atmospheric composition processes are discussed for US (Zhang, 2008) and European (Baklanov et al., 2014) models.Klein et al. (2012) extended applications of coupled models for "biological weather", defined as "the short-term state and variation of concentrations of bioaerosols", in particular for pollen modelling and forecasting.
The online integration of meso-meteorological models (MetM) and atmospheric aerosols and ACT models gives a possibility to utilise all meteorological 3-D fields in the ACT model at each time step and to consider nonlinear feedbacks of air pollution (e.g.atmospheric aerosols) on meteorological processes/climate forcing and then on the chemical composition of the atmosphere.This very promising way for future atmospheric modelling systems (as a part of, and a step toward, the Earth system modelling, ESM) will lead to a new generation of seamless coupled models for meteorological, chemical and biochemical weather forecasting.A seamless approach for "one atmosphere" integrated meteorology-chemistry/aerosol forecasting systems is analysed by the COST Action ES1004 EuMetChem (see, e.g.Baklanov et al., 2015), and an overview of the current state of online coupled chemistry-meteorology models and needs for further developments were published in Zhang (2008), WMO (2016), Baklanov et al. (2017) and Sokhi et al. (2017).
The methodology on how to realise the suggested integrated concept was demonstrated on a European example of the Enviro-HIRLAM (Environment -High Resolution Limited Area Model) integrated modelling system (Baklanov et al., 2008a;Korsholm, 2009).Experience from first HIRLAM community attempts to include pollutants into the NWP model (Ekman, 2000) and from pioneering online coupled meteorology-pollution model developments of the Novosibirsk science school (Marchuk, 1986;Penenko and Aloyan, 1985;Baklanov, 1988) was actively used for developments of the Enviro-HIRLAM modelling system.
The Enviro-HIRLAM is developed as a fully online integrated NWP and ACT modelling system for research and forecasting of meteorological, chemical and biological weather.The integrated modelling system is developed by the Danish Meteorological Institute (DMI) and other collaborators (Chenevez et al., 2004;Baklanov et al., 2008aBaklanov et al., , 2011b;;Korsholm et al., 2008Korsholm et al., , 2009;;Korsholm, 2009) and included as the baseline system of the chemical branch of the HIRLAM consortium (Fig. 1).
The model development was initiated at DMI more than 15 years ago and it is used now in several countries.The modelling system is being used for different completed and ongoing research projects (FP6 FUMAPEX; FP7 MEGAPOLI, PEGASOS, MACC, TRANSPHORM, MarcoPolo; Nord-Forsk NetFAM, MUSCATEN, CarboNord, CRAICC-PEEX, CRUCIAL; COST Actions -728, 732, ES0602 ENCWF, ES1004 EuMetChem) and has been used for operational pollen forecasting in Denmark since 2009 (Rasmussen et al., 2006;Mahura et al., 2006b) and operational atmospheric composition (with a focus on aerosols) for China since November 2016 (Mahura et al., 2016(Mahura et al., , 2017)).Following the main strategic plans (Baklanov, 2008;Baklanov et al., 2011a) within HIRLAM-B and -C projects, further developments of the modelling system will be shifting to a new NWP platform (from HIRLAM to HARMONIE), and a close collaboration with the ALADIN (Aire Limitée Adaptation dynamique Développement InterNational) community was initiated in 2014.
In this paper, an overall description of the current version of the Enviro-HIRLAM coupled modelling system with improved parameterisations of meteorology-composition twoway interactions, the main steps in its development and examples in different application areas for air quality, weather and pollen forecasting are considered for the first time.Section 2 provides a detailed description of the Enviro-HIRLAM modelling system and its key developments in the meteorological core, chemistry and aerosol dynamics parts, aerosol-meteorology interactions, models urbanisation and improvements of numerical algorithms.Section 3 describes a few types of Enviro-HIRLAM applications for meteorological and environmental forecasting and assessment studies.Sections 4 and 5 continue discussions and summarise the model applicability and provide recommendations for future research.Appendix A includes brief information about the Enviro-HIRLAM model development history.A list of acronyms is provided in Appendix B.

Modelling system structure
The Enviro-HIRLAM is a fully online coupled (integrated) NWP and ACT modelling system for research and forecasting of meteorological, chemical and biological weather (see schematics in Fig. 2).The modelling system was originally developed by DMI and further with other collaborators, and now it is included by the European HIRLAM consortium as a baseline system in the HIRLAM Chemical Branch (http: //hirlam.org/index.php/documentation/chemistry-branch).It was the first mesoscale online coupled model in Europe that considered two-way indirect feedbacks between meteorology and chemistry/aerosols (WMO-COST, 2008).
The main steps of the model development were realised, including (i) model nesting for high resolutions, (ii) improved resolving planetary boundary layer (PBL) and surface layer structure, (iii) urbanisation of the NWP model, (iv) improvement of advection schemes, (v) emission inven-tories and models, (vi) implementation of gas-phase chemistry mechanisms, (vii) implementation of aerosol dynamics, (viii) realisation of aerosol feedback mechanisms.
The first version was based on the DMI-HIRLAM NWP model with online integrated pollutant transport and dispersion (Chenevez et al., 2004), chemistry, deposition and indirect effects (Korsholm et al., 2008;Korsholm, 2009) and later aerosol (only for sulfur particles) dynamics (Baklanov, 2003;Gross and Baklanov, 2004).To make the model suitable for chemical weather forecasting in urban areas, the meteorological part was improved by implementation of urban sublayer parameterisations (Baklanov et al., 2008b;Mahura et al., 2008a;González-Aparicio et al., 2013).The model's dynamic core was improved by adding a locally massconserving semi-Lagrangian numerical advection scheme (Kaas, 2008;Sørensen, 2012;Sørensen et al., 2013), which improves forecast accuracy and enables performing longer runs.More details of the system development history are presented in Appendix A.
The current version of Enviro-HIRLAM (Nuterman et al., 2013(Nuterman et al., , 2015) ) is based on the Reference-HIRLAM v7.2 with a more advanced and effective chemistry scheme, multi-compound modal approach aerosol dynamics modules, aerosol feedbacks on radiation (direct and semi-direct effects) and on cloud microphysics (first and second indirect effects).This version is continuously under development and evaluation for various weather and air-quality-related applications (in particular, within the COST Action ES1004 where the above-mentioned effects were extensively discussed; see, e.g.Baklanov et al., 2014).
Vertical and horizontal resolutions of the model are flexible.Limitations, e.g.due to the hydrostatic approximation, are provided (minimum 1.5 km of the horizontal resolution for flat terrains, e.g. for Copenhagen).

Meteorological core of the system
The first version of Enviro-HIRLAM was based on a previous version HIRLAM-tracer.At its meteorological core lies DMI-HIRLAM, version 6.3.7,employed for limited-area short-range operational weather forecasting at DMI (Chenevez et al., 2004).The current model version used in studies is based on the reference version of the HIRLAM community meteorological NWP model HIRLAM version 7.2 and an online coupled environmental block (called the Enviro-) allowing to take into account spatial-temporal evolution of atmospheric chemical and biological aerosols driven by meteorology from the NWP block.
HIRLAM is a hydrostatic NWP model which is used for both research and operational purposes.The model provides a forecast of the main meteorological fields, including air temperature and specific humidity, atmospheric pressure, wind speed and direction, cloud cover and turbulent kinetic energy (TKE) based on forward-in-time integration of the primitive equations (dynamical core) (Holton, 2004) and physical processes such as radiation, vertical diffusion, convection, condensation, etc. (physical core).
The detailed NWP HIRLAM description can be found in the HIRLAM reference guide scientific documentation (Undén et al., 2002) and its following upgrades and modifications (for more details, see http://www.hirlam.org).
The hydrostatic approximation of the model can be a limitation for increasing the resolution for urban simulations.However, sensitivity tests for a medium-size city demonstrated that 2.5 km was the optimal resolution, allowing at the same time to obtain satisfactory reproducibility of the largescale processes and to explore the urban effects at local scale (González-Aparicio et al., 2013).For other metropolitan areas, such as Paris, Rotterdam, St. Petersburg and Shanghai, a similar resolution was chosen, whereas for Copenhagen (with its flat terrain) the highest suitable resolution tested was 1.5 km and provided reasonable verification results (Mahura et al., 2006a(Mahura et al., , 2008b(Mahura et al., , c, 2016)).Within a selected metropolitan area, there could be only a few grid cells having 100 % representation of the urban fraction, but taking into account all urban grid cells, the boundaries of the cities (number of cells) could be substantially larger.Moreover, most of the existing parameterisations in the physics core of any NWP model might need a revision when resolutions of 1 km and finer are used.
Following the main strategic development within HIRLAM (HIRLAM-B and -C projects), there are plans for further developments of Enviro-HIRLAM shifting to a new non-hydrostatic NWP platform (e.g.HARMONIE model) and incorporating chemistry modules and aerosolradiation-cloud interactions into the future integrated system (Baklanov, 2008;Baklanov et al., 2011a).
The new non-hydrostatic version under HARMONIE is under development and only some elements are realised so far.The non-hydrostatic HARMONIE-AROME model includes only some aerosol effects.The physics included in this version of HARMONIE has recently been detailed by Bengtsson et al. (2017).HARMONIE-AROME is based partly on Meso-NH (mesoscale non-hydrostatic atmospheric model), which is a cloud-resolving model that includes stateof-the-art chemistry and aerosol interactions (e.g.Berger et al., 2016).However, Meso-NH cannot be run as a near-realtime NWP model, as is possible with Enviro-HIRLAM.

Tropospheric sulfur cycle
The simple tropospheric sulfur cycle chemistry module in Enviro-HIRLAM, used for long-term runs (up to 1 year), is based on the sulfur cycle mechanism developed by Feichter et al. (1996) treating three prognostic species: dimethyl sulfide (DMS), sulfur dioxide (SO 2 ) and sulfate (SO 2− 4 ).The mechanism includes DMS and SO 2 oxidation by hydroxyl (OH) and DMS reactions with nitrate radicals (NO 3 ) in the gas-phase part.The heterogeneous aqueous phase chemistry is comprised of SO 2 oxidation reactions by H 2 O 2 and O 3 .Accounting for dissolution effects of SO 2 in the aqueous phase is performed according to Henry's law.An output of the global chemistry transport model MOZART (Horowitz et al., 2003) is used to prescribe three-dimensional oxidant fields of OH, H 2 O 2 , NO 2 and O 3 .
The sulfate produced in the gas phase is referred to the gases and can be condensed on pre-existing aerosols or to nucleate by the aerosol microphysics M7 module (see Sect. 2.4).Moreover, in-cloud produced sulfate is accumulated on the pre-existing accumulation and coarse-mode aerosols.
The tropospheric sulfur cycle chemistry is used together with the M7 aerosol microphysics module because of its relative simplicity and low computational cost.The Carbon Bond Mechanism version Z (CBM-Z) gas-phase chemistry (see the next section) is not interfaced with the M7 aerosol module because of several reasons: (1) the aerosol microphysics module does not include secondary organic aerosols; therefore, there is no need for a complex gas-phase mechanism with volatile-organic-compound-related reactions, and (2) it is too computationally expensive to use CBM-Z together with M7 for both weather and atmospheric composition prediction.

Gas-phase chemistry
The gas-phase chemistry scheme consists of sets of chemical schemes running from simple schemes for chemical weather forecasts to highly complex schemes for research and case studies.In order to make the model computationally efficient for different applications and operational forecasting, several condensed atmospheric chemical schemes have been tested into Enviro-HIRLAM since the first version of the model system was realised (Korsholm, 2009;Gross and Baklanov, 2004).In the current version of Enviro-HIRLAM, the tropospheric condensed CBM-Z (Zaveri and Peters, 1999), a variant of the CBM-IV gas-phase chemistry scheme (Gery et al., 1989), with a fast solver based on radical balances (Sandu and Sander, 2006), has been implemented in the model.CBM-Z uses lumped species that represent broad categories of organics based on carbon bond structure.It is closely related to CBM-IV which is widely used in air pollution eval-uations but with expansions to include reactions that are important in the remote troposphere.It also uses the most general organic category (PAR for paraffin) to represent miscellaneous carbon content so that carbon mass is conserved.
Six environmental/smog chamber experiments were used to validate the gas-phase schemes as box models and within a regional climate model (Shalaby, 2012;Shalaby et al., 2012).The Tennessee Valley Authority (TVA) and the EPA chamber experiments were used to evaluate the different gas-phase schemes and different chemical solvers.Namely, TVA005 and TVA006 are designed to test the simple system of NO x ; TVA068 is designed to test a simple mixture of volatile organic compounds (VOCs) with very high NO x .EPA069A, EPA073A and EPA150A are used to validate the schemes with low NO x concentration and high VOC concentration.

Chemical solvers
Calculating the time evolution of gas-phase chemistry requires a numerical integration of a set of stiff ordinary differential equations (ODEs) and is among the most computationally expensive operations performed in a photochemical grid model.The equations for photochemical production and loss are computationally expensive because they form a stiff numerical system.The photochemical mechanisms described above were implemented using two different chemical solvers to solve the tendency equation for photochemical production and loss: (1) the Rosenbrock (ROS) solver (Sandu et al., 1997;Hairer and Wanner, 1996) as implemented by the Kinetic PreProcessor (KPP) (Sandu and Sander, 2006) and (2) the computationally rapid radical balance method (RBM) of Sillman (1991).RBM utilises the fact that much of the complexity of tropospheric chemistry stems from the HO x radical family (OH, HO 2 and RO 2 ), which has a limited set of sources and sinks.The method solves reverse Euler equations for OH and HO 2 based on the balance between sources, sinks and (if applicable) prior concentrations at the start of the time step.Reverse Euler equations for other species are solved in a reactant-to-product order, in some cases involving pairs of rapidly interacting species, and with some modifications to increase accuracy in exponential decay situations.The procedure is equivalent to a reverse Euler solution using sparse-matrix techniques but with the matrix inversion linked specifically to the behaviour of OH and other species in the troposphere.Prior work tested several atmospheric chemistry mechanisms in the model by taking into account different chemical solvers.We select the photochemical mechanism CBM-Z because it affords a reasonable trade-off between accuracy and computational efficiency.During the prior work, including the validation stages of the gas-phase schemes (results not shown), we used KPP to generate the Fortran code of three different gas-phase schemes: CBM-Z (Zaveri and Peters, 1999), GEOS-CHEM (Evans et al., 2003) and the Regional Atmospheric Chemistry Model (RACM; Stockwell et al., 1997).In order to fit within our main aim of the chemi-cal weather predication, we did not use both GEOS-CHEM and RACM because they are very computationally expensive schemes due to their extensive number of chemical reactions.The KPP provides a flexible tool to generate a wellcoded chemical mechanism according to the user choice of a given ODE solver.We use KPP tools to create the gas-phase chemical mechanisms including the solvers for three chemical mechanisms.Usually, the Rosenbrock solver is selected for most of simulations due to its ability to be a fast computational solver (Sandu et al., 1997).

Photolysis rates
Photolysis rates are determined as a function of various meteorological and conditional inputs.Rates for specific conditions are determined by interpolating from an array of predetermined values.The latter is based on the Tropospheric Ultraviolet-Visible (TUV) model developed by Madronich and Flocke (1999), using a pseudo-spherical discrete ordinate method (Stamnes et al., 1988) with eight streams.The eight-stream TUV is the most accurate method for determining photolysis rates but is computationally too expensive for use in 3-D models.Photolysis rate constants are calculated using the Fast-J radiative transfer model (Wild et al., 2000) with O( 1 D) quantum yields updated to JPL2003 (Sander et al., 2003).
For simplicity, photolysis rates are estimated as the following.At first, for the simple reactions, the photolysis rates are estimated as a function of the number of parameters such as meteorological and chemical inputs including altitude, solar zenith angle, overhead column densities for O 3 , SO 2 and NO 2 , surface albedo, aerosol optical depth, aerosol single scattering albedo, cloud optical depth and cloud altitude.Then, for the complex reactions, the photolysis rates are estimated as a lookup table using the TUV model.TUV is run offline and used to calculate a lookup table of the photolysis rates, and then this lookup table is implemented under different weather conditions inside the model.
Photolysis rates can be significantly affected by the presence of clouds.Cloud optical depths are determined using the random overlap treatment described by Feng et al. (2004), which assumes that cloudy and cloud-free subregions in each model grid box randomly overlap with cloudy and cloud-free subregions in grid boxes located above or below (Briegleb, 1992).The method used to correct for cloud cover is based on Chang et al. (1987), which requires information on cloud optical depth for each model grid cell.Optical depth is used to reduce photolysis rates for layers within or below clouds to account for UV attenuation or to increase photolysis rates due to above-cloud scattering.In general, below-cloud photolysis rates will be lower than the clear-sky value due to the reduced transmission of radiation through the cloud.Similarly, photolysis rates are enhanced above the cloud due to the reflected radiation from the cloud.Cloud optical depths and cloud altitudes from Enviro-HIRLAM are used in the pho-tolysis calculations, thereby directly coupling the photolysis rates and chemical reactions to meteorological conditions at each model time step.

Heterogeneous chemistry
Many gas-phase species are water soluble, and sulfate and ammonia together with water take part in binary/ternary nucleation.In order to consider these processes, a simplified liquid-phase equilibrium mechanism with the most basic equilibria is included in NWP-Chem-Liquid.The NWP-Chem-Liquid is a thermodynamic equilibrium model, described in Korsholm et al. (2008).This equilibrium module is solved using the analytical equilibrium iteration method (Jacobson, 1999).The reactions are summarised in Korsholm (2009) and the module will be updated to include the impact of organic compounds from anthropogenic and biogenic sources.

Aerosol dynamics module
The first aerosol module in Enviro-HIRLAM was based on the CAC (chemistry-aerosol-cloud) model with the modal approach for description of aerosol size distribution (Baklanov, 2003;Gross and Baklanov, 2004) and considered only sulfur-type aerosols (Korsholm, 2009).
The current version of the Enviro-HIRLAM model has an M7 aerosol microphysics module (Vignati et al., 2004) together with aerosol removal processes ported from the ECHAM5-HAM climate model (Stier et al., 2005).There are two types of particles considered: insoluble and mixed (water-soluble) particles.The particles are split into seven classes depending on particle size and solubility by means of the "pseudo-modal" approach.Four classes are used to represent mixed particles, i.e. nucleation, Aitken, accumulation and coarse modes, and another three classes are for the insoluble (Aitken, accumulation and coarse modes).Four predominant aerosol types are included -black carbon (BC) and primary organic carbon (OC), sulfate, mineral dust and sea salt.The M7 aerosol dynamics includes nucleation, coagulation and sulfuric acid condensation processes.Coagulation and condensation will lead to formation of mixed particles from the insoluble ones.Different aerosol types mentioned above (as well as others, e.g.pollen particles) are provided as separate species in the model outputs along with lumped PM 10 and PM 2.5 .

Dry deposition and sedimentation
The dry deposition fluxes of gases and aerosols (for both number and mass concentrations) are calculated from the aerodynamic, quasi-laminar boundary layer as the product of the surface layer concentration and the dry deposition velocity (Stier et al., 2005).The fluxes are used as the lower boundary condition in the semi-implicit vertical diffusion TKE-CBR scheme (Cuxart et al., 2000).The calculation of the dry deposition velocities is performed by means of the serial resistance approach.The "big-leaf" method is used to calculate surface resistance (Ganzeveld and Lelieveld, 1995;Ganzeveld et al., 1998) per each grid cell for the snow/ice, water, bare soil, low-vegetation and forest surface types.The SO 2 soil resistance is a function of soil pH, relative humidity, surface temperature and the canopy resistance, while surface resistances for other gases are prescribed.The canopy resistance is computed from stomatal resistance and monthly mean leaf area index (LAI) values from the Enviro-HIRLAM interaction soil-biosphere-atmosphere scheme (Noilhan and Planton, 1989).
The sedimentation of the aerosol particles is calculated throughout the atmospheric column.The calculation of the sedimentation velocity is based on the Stokes velocity with the Cunningham slip-flow correction factor accounting for non-continuum effects (Seinfeld and Pandis, 2006).In order to satisfy the Courant-Friedrichs-Lewy stability criterion, the sedimentation velocity is limited by the ratio of the model layer thickness and the time step.

Wet deposition
There are several options for the wet deposition in the model.The first version used the aerosol-size-dependent parameterisation of Baklanov and Sørensen (2001).In the latest version, fixed size-and composition-dependent scavenging parameters are also applied for wet deposition calculation and are different for stratiform and convective clouds (Stier et al., 2005).They were derived from measurements of interstitial and in-cloud aerosol contents.These scavenging coefficients depend on the aerosol modes, total cloud water and fraction (liquid and ice) and the conversion rates of cloud liquid water and cloud ice to precipitation through autoconversion, aggregation and accretion processes.The precipitation re-evaporation before it reaches the ground is also included.The Soft TRAnsition COndensation (STRACO) cloud scheme (Sass, 2002) provides water and ice precipitation fluxes, normalised by the precipitation rates, to the wetdeposition scheme, which uses prescribed size-dependent collection efficiencies for rain and snow (Seinfeld and Pandis, 1998).

Emission modules and preprocessor
The model includes anthropogenic, biomass burning (wildfires) and natural emission fluxes of both gases and aerosols.These emissions are processed in different ways, because some of them are pure datasets derived from ground-based and satellite observations.The others are interactively developing during the model integration and depend on the meteorological conditions at the current time step and land-use, land-cover or water surface types.The anthropogenic emis-sion inventory developed by the Netherlands Organisation for Applied Scientific Research (TNO) (Kuenen et al., 2014).Linked to the model is a dataset of yearly accumulated fluxes of gases, such as CO, CH 4 , NO x , SO 2 , NH 3 , non-methane volatile organic compounds (NMVOCs) and particulate matter (PM) in two size bins (2.5 and 10 µm) which are attributed to 10 source sectors, including energy industries, residential combustion, industry, etc., denoted by SNAP (Selected Nomenclature for sources of Air Pollution) codes.The inventory has a resolution of 0.06 • × 0.12 • and covers the entirety of Europe, the European part of Russia, north of the Sahara and a part of the Middle East.Total NMVOC emissions are split into 25 VOC compound groups by source sectors by country (Kuenen et al., 2010).The PM 2.5 and PM 10 emissions, split into six aerosol species (BC, OC, Na, SO 4 , other coarse primary and other fine primary particles), are applied following TNO recommendation (Kuenen et al., 2010).Because the dataset contains accumulated surface fluxes, one needs to redistribute them in order to reproduce diurnal, weekly and monthly emissions variability.The emissions can also occur at different heights; e.g.emissions from power plants are elevated and those from traffic are at the surface, and thus vertical redistribution is applied within the first eight model hybrid levels.Therefore, temporal and vertical profiles developed by TNO for different gaseous and aerosol species and SNAP codes are used in the emission preprocessor.The global biomass burning (wildfire) emission inventory (IS4FIRES; Sofiev et al., 2012), developed by the Finnish Meteorological Institute (FMI), has similar structure except for the number and kinds of available gaseous and aerosol species as well as the resolution.The inventory data are the total PM flux.The flux is split into PM 2.5 and coarse PM consisting of ash.The PM 2.5 primarily consists of organic and black carbon (OC and BC) and a remainder of organic matter that is not carbon; for details, see Andreae and Merlet (2001).The biomass burning emissions typically show a diurnal cycle variability, and therefore corresponding coefficients are applied (Giglio, 2007).The wildfire emissions are also redistributed vertically with different proportions -the lowest at 200 m and the highest up to 1 km over the ground.
The natural emissions of gases and aerosols are fully interactive and calculated online.There is dimethyl sulfide (DMS; Nightingale et al., 2000) emission from oceans, which depends on the wind speed and seasonal variability of DMS solution in the water.Soluble sea-salt aerosol emissions (Zakey et al., 2008) are driven by wind speed and temperature, and insoluble mineral dust aerosol emissions (Zakey et al., 2006) also depend on meteorology as well as hydrological parameters.Both sea-salt and dust aerosols are emitted in accumulation and coarse modes.

Direct and semi-direct effects
Enviro-HIRLAM contains parameterisations of the direct and semi-direct effects of aerosols.Direct and semi-direct effects are realised by modification of the Savijärvi radiation scheme (Savijärvi, 1990;Wyser et al., 1999) with implementation of new fast analytical short-wave (SW) and long-wave (LW) aerosol transmittances, reflectances and absorptances.The two-stream approximation equations for anisotropic non-conservative scattering described by Thomas and Stamnes (2002) are used for these calculations.The Global Aerosol Data Set/Optical Properties of Aerosols and Clouds (GADS/OPAC) aerosols of Köpke et al. (1997) are used as input to the routine.The species include BC (soot), minerals (nucleus, accumulation, coarse and transported modes), sulfuric acid, sea salt (accumulation and coarse modes), "water soluble" and "water insoluble" aerosols.In addition to the more standard nucleation, accumulation and coarse aerosol size modes, we consider, according to Köpke et al. (1997), the transported size mode to describe aerosols that have been transported over a long distance, for instance, Saharan aerosols that have been blown to the Atlantic Ocean.In order to make the calculations fast, optical properties that are spectrally averaged over the entire SW and LW spectra are used.The spectra used are shown in Fig. 3.The shortwave spectrum is a clear-sky spectrum from 2 km height in a standard atmosphere (Anderson et al., 1986) calculated with the DISORT algorithm (Stamnes et al., 1988) run in the Li-bRadtran framework (Mayer and Kylling, 2005).The longwave spectrum is calculated similarly and is based on the overall atmospheric LW transmittance of a standard atmosphere.

First and second indirect effects
For aerosol-cloud interactions, a modified version of the STRACO cloud scheme (Sass, 2002) is used in Enviro-HIRLAM.This scheme developed for operational NWP has recently been upgraded using new efficient methods to account for aerosol effects on cloud formation and microphysics.The scheme is able to account for convective transports of new variables.The prognostic aerosol fields are coupled directly to the cloud physical and microphysical properties.Liquid cloud droplet number is calculated based on aerosol size, number and solubility, and the STRACO subgrid supersaturation field is used as the basis for the droplet nucleation calculation.This ensures consistency with the cloud water mass.
The modelled liquid droplet number evolves in time according to the following processes: droplet nucleation, selfcollection, sedimentation and evaporation.In order to close the tendency calculations, the liquid cloud droplet distribution is assumed to follow a gamma distribution where  (2000).Self-collection is the process whereby droplets collide and stick together but do not become raindrops.The parameterisation of self-collection processes follows Seifert and Beheng (2006).Sedimentation is calculated to be consistent with the mass of rainwater in a given model time step under the basic assumption that the largest droplets are removed first from the cloud.Similarly, evaporation of a droplet below activation radius is calculated to be consistent with the total evaporated cloud water under the assumption that the smallest droplets evaporate first.
The cloud droplet effective radius controls the liquid phase absorptivity and transmissivity and is calculated from liquid water mass and droplet number.Here, it is also dependent on the shape of the droplet distribution which evolves in time.Autoconversion follows Rasch and Kristjansson (1998) and is directly dependent on the calculated droplet number.
Abdul-Razzak and Ghan (2000) parameterisation for aerosol activation has been extensively tested in many online coupled weather and climate models.However, the STRACO cloud microphysics scheme with parameterisations of aerosol activation, cloud droplet nucleation, sedimentation, evaporation and self-collection has been evaluated only with 1-D column HIRLAM, so it needs to be further thoroughly evaluated.

Urban parameterisations and models urbanisation
The representation of urban areas in Enviro-HIRLAM contains the following aspects and processes (Baklanov et al., 2005): i. model down-scaling, including increasing vertical and horizontal resolution and nesting techniques; ii. modified high-resolution urban land-use classifications, parameterisations and algorithms for roughness parameters in urban areas based on the morphologic method; iii. specific parameterisation of the urban fluxes in the mesoscale model; iv.modelling/parameterisation of meteorological fields in the urban sublayer; and v. calculation of the urban mixing height based on prognostic approaches.
The nesting technics and downscaling methods are actively and successfully used for urban areas to reach the necessary resolution for resolving or parameterisation of urban features and effects.The details of this approach with the Enviro-HIRLAM model were described, e.g. in Baklanov and Nuterman (2009).With respect to metropolitan areas, the downscaling for finer resolution allows to reproduce smaller-scale meteorological patterns, and then these patterns are further modified through running urban parameterisation modules only for grid cells where the cities are presented.
The urban parameterisations in the model contain three different approaches which may be combined.The first is the simplest implementation which contains modifications of the surface roughness, the anthropogenic heat flux, the storage heat flux and the albedo over urban areas.These are identified in the model using urban fractions extracted from the land-use database (CORINE) employed at DMI (Mahura et al., 2005b(Mahura et al., , 2006a(Mahura et al., , 2007a;;Baklanov et al., 2005Baklanov et al., , 2008b)).The first module is the computationally cheapest way of urbanising the model and it can be used for operational NWP as well as for regional climate modelling.The second modulebuilding effect parameterisation (BEP) (Martilli et al., 2002) -gives a possibility to consider the energy budget components and fluxes inside the urban canopy although it is relatively more expensive (5-10 % computational time increase) (Mahura et al., 2008b(Mahura et al., , c, 2010b;;Fig. 4).However, this approach is sensitive to the vertical resolution of NWP modwww.geosci-model-dev.net/10/2971/2017/Geosci.Model Dev., 10, 2971-2999, 2017 2980 A. Baklanov et al.: Enviro-HIRLAM online integrated meteorology-chemistry modelling system els and is not very effective if the first model level is higher than 30 m.Therefore, the increasing of the vertical resolution of current NWP models is required.The third module -Soil Model for Submesoscales, Urbanized (SM2-U) version (Dupont and Mestayer, 2006;Dupont et al., 2006) -is considerably more expensive computationally than the first two modules (Mahura et al., 2005a;Baklanov et al., 2008b).
However, the third one provides the possibility to accurately study the urban soil and canopy energy exchange including the water budget.Therefore, the BEP scheme is considered as the baseline option and third SM2-U module is recommended only for use in advanced urban-scale NWP and meso-meteorological research models.The details of implementations of different urban modules, own developments and comparisons of different approaches and modules were published in previous papers (Mahura et al., 2005a(Mahura et al., , b, 2006a(Mahura et al., , 2008a(Mahura et al., , b, c, 2010b;;Baklanov et al., 2005Baklanov et al., , 2008b)).The main approach includes an integration of the urban modules into the ISBA (interaction soil-biosphere-atmosphere) land surface scheme of the NWP/HIRLAM model.The urban modules are activated only on those grid cells of the model domain where the urban fraction is presented.
The urban boundary layer is very inhomogeneous and plays an important role in forming urban meteorological fields and especially in dispersion of atmospheric pollutants.Therefore, for calculation of the urban mixing height, in addition to the common diagnostic approaches, prognostic equations were used according to Zilitinkevich et al. (2002) and Zilitinkevich and Baklanov (2002).

Transport schemes
Until 2012, there were basically two options for transport schemes in Enviro-HIRLAM (Chenevez et al., 2004): (a) the traditional non-conserving but highly efficient semi-Lagrangian (SL) scheme (Robert, 1981) in HIRLAM and (b) the much less efficient flux-based and positive definite finite volume scheme by Bott (1989) with updates by Easter (1993).In 2012, the default transport scheme was updated to a new monotonic version of the locally massconserving semi-Lagrangian (LMCSL) scheme (Kaas, 2008;Sørensen et al., 2013).This scheme, used in the present version of Enviro-HIRLAM, to be described briefly below, is almost as efficient as the traditional SL scheme but now has the attractive properties of inherent mass conservation, and it is monotonic and positive definite.
In HIRLAM and former versions of Enviro-HIRLAM, a traditional SL scheme is used for advecting the specific concentration of water constituents or the mixing ratio q i of any tracer i. Considering mixing ratio, this means that when ignoring any sources/sinks and turbulent mixing the prognostic transport equation to be solved is simply The traditional SL numerical integration of Eq. ( 1) reads where subscript k is the grid point/cell index and superscripts n and n + 1 represent two consecutive time steps, respectively.The subscript * k indicates the tri-cubic interpolation to the location of the departure point of the upstream trajectory, which arrives in grid point k at time level n + 1.The tri-cubic interpolation in Eq. ( 2) can also be represented as a sum of interpolation weights involving 64 grid points surrounding the departure point.Formally, this can be expressed as where K is the total number of grid points in the entire integration domain.Note that for each k only 64 w k,l weights are different from zero.When converting mixing ratio into volume density, i.e. (ρ i ) n+1 k = (ρ d ) n+1 k (q i ) n+1 k , and subsequently summing over the integration area, the traditional SL scheme is not mass conserving.Therefore, in LMCSL (Kaas, 2008), a different approach is followed, namely, as in most other mass-conserving transport schemes, to solve the complete continuity equation: still omitting sources/sinks and turbulent mixing and then evaluating the mixing ratio from In LMCSL, Eq. ( 4) is solved in a rather unusual way by modifying the interpolation weights in Eq. ( 3) in such a way that the sum of mass given off at time step n by a Eulerian grid cell l to all departure points that it influences is exactly equal to its own mass.In other words, LMCSL is based on simple partition of unity.The modified weights become where V k is the volume of Eulerian grid cell k.Using the modified weights the basic LMCSL forecast reads as follows: As the traditional SL scheme, the LMCSL is not inherently monotonic or positive definite.Therefore, an a posteriori iterative locally mass-conserving (ILMC) filter was developed (Sørensen et al., 2013).This filter ensures that the mixing ratio of the forecast will never be larger/smaller than the largest/smallest mixing ratio of the eight grid cells surrounding the upstream trajectory departure point at time level n.The ILMC filter designed to be as local as possible since nonlocal filters will generate non-physical chemical reactions.This is ensured by an iterative approach where the mass discrepancy is redistributed among the neighbouring cells in the first iteration, and the distribution radius is increased, in the case that there is remaining mass discrepancy, for the next iteration(s).In general, one or two iteration(s) are sufficient.
The LMCSL transport scheme in combination with the ILMC produces accurate monotonic and positive definite forecasts for water vapour, liquid/ice water and chemical constituents.As an example, the simulated PM 2.5 concentration on 17 July 2010 with horizontal resolution of approximately 16 km is shown in Fig. 5.It can be seen that the model is able to reproduce, e.g.sharp transitions related to fronts over the North Atlantic.A more in-depth analysis of the ability of ILMC to reproduce sharp gradients can be found in Sørensen et al. (2013), in particular Fig. 3 and the accompanying discussion in that paper.
It should be noted that the dynamical core in Enviro-HIRLAM is identical to that of HIRLAM.Thus, the dry-air density for dynamics is calculated using a traditional SL approximation to Eq. ( 4), i.e. not the LMCSL.Therefore, the Enviro-HIRLAM is not formally wind-mass consistent regarding tracer transport.However, the large-scale precipita-tion fields in the traditional HIRLAM and Enviro-HIRLAM are very similar (see, e.g.Fig. 4 in Sørensen et al., 2013), which suggests that wind-mass inconsistency is of minor importance.In principle, no monotonic transport schemes can be wind-mass consistent since the monotonic limiters formally destroy the consistency (see discussions on the issue of wind-mass inconsistency in atmospheric models in Jöckel et al., 2001).

Modelling system applications
Possible applications of the online integrated Enviro-HIRLAM modelling system include the following: chemical weather forecasting, air quality and chemical composition longer-term assessment, weather forecast (e.g. in urban areas, severe weather events), pollen and bio-aerosol transport forecasting, climate change modelling, studies of climate change effects on atmospheric pollution on different scales, anthropogenic impacts on atmospheric processes, weather modifications, geo-engineering, contamination from volcano eruptions, sand and dust storms, nuclear explosion consequences and other emergency preparedness modelling.
Several realised/tested types of applications of the Enviro-HIRLAM for meteorological, environmental and climate forecasting and assessment studies are highlighted in Fig. 1 and will be demonstrated below.

Applications for numerical weather prediction
Several Enviro-HIRLAM sensitivity and validation studies of aerosol feedbacks on meteorological processes were done previously (see, e.g.Korsholm, 2009;Korsholm et al., 2010;Baklanov et al., 2011a, b;Sokhi et al., 2016).For example, the effects of urban aerosols on the urban boundary layer height, can be comparable to the effects of the urban heat island ( h is up to 100-200 m for stable boundary layer) (Baklanov et al., 2008a).Further studies (Korsholm et al., 2010) of megacity effects on the meteorology/climate and atmospheric composition showed that aerosol feedbacks through the first and second indirect effects induce considerable changes in meteorological fields and large changes in chemical composition (see Sect. 3.4) in the case of convective clouds and little precipitation.The monthly averaged changes in surface temperature due to aerosol indirect effects of primary aerosol emissions in western Europe were analysed and validated versus measurement data.It was found that a monthly averaged signal (difference between runs with and without the indirect effects) in surface temperature can reach 0.5 • C (Fig. 2.2b in Korsholm et al., 2010).Korsholm (2009) studied the impact of aerosol indirect effects on surface temperatures and air pollutant concentrations for a 24 h simulation over a domain in northern France, including Paris, in a convective case with low precipitation.He found a marginally improved agreement with observed 2 m temperatures and a marked redistribution of NO 2 in the domain, primarily as a result of the second indirect effect.
To perform analysis of atmospheric aerosol effects on clouds and precipitation, the year 2010 was selected for Enviro-HIRLAM simulations.That year, especially summer, was characterised by severe weather events such as floods, heat waves and droughts across the Middle East, most of Europe and European Russia.The model was forced by boundary and initial conditions produced by the ECMWF IFS (IFS-CY40r1) and MOZART (Horowitz et al., 2003) models for meteorology and atmospheric composition, respectively.The Enviro-HIRLAM modelling domain with horizontal resolution of 0.15 • × 0.15 • , 310 × 310 grid cells and 40 vertical hybrid sigma levels extending to pressures less than 10 hPa covers Europe, north of the Sahara and European Russia.The modelling domain was partitioned into 120 CPU cores and the model was run with a time step of 300 s.The model includes emissions from anthropogenic sources developed by TNO and from wildfires produced by FMI as well as interactive DMS, sea salt and dust emissions (for details, see Sect.2.5).
For aerosol-cloud interactions, these were estimated also for July 2010 by means of delta functions, i.e. the difference between outputs of models: Enviro-HIRLAM with aerosolcloud interactions (ENV) and Reference-HIRLAM (REF).Figure 6a shows deltas (ENV-REF) of total cloud cover over model domain, which is mainly increased (with local maxima up to 90 %) except several inland areas, such as Finland, the borders of Germany, Poland and Austria, where cloud cover decreased by almost 10-fold.The ENV runs revealed the increase of average cloud top height by approximately 2 %.The delta function of cloud water content at average cloud base shows (Fig. 6b) its increase compared to REF and local maxima over the North Atlantic, North Sea, Sweden, Switzerland and Austria.These areas are occupied by precipitating clouds as seen in Fig. 7.
The absolute frequencies of stratiform and convective precipitation over computational domain are decreased compared to the REF model, while the amount of convective precipitation during heavy precipitation events is increased.Hence, the wet deposition of particles decreases in summer because it rather depends on precipitation frequency than on its amount.The REF model run tends to overpredict both frequency and amount of precipitation.But the inclusion of aerosol-cloud interactions can improve general model performance; i.e. the ENV run bias for precipitation with respect to its frequency and amount has been decreased compared to the REF model run (Fig. 8).
Sensitivity studies on the model response to aerosol effects indicate strong "signals", but they do not guaranty improvements.For example, Korsholm (2009) considered evaluations only for some elements (e.g. the coupling interval) in the previous analysis and made corresponding conclusions about the improvements.Other feedback mechanisms, especially for aerosol-cloud interactions, were analysed mostly as sensitivity studies or evaluations for short-term episodes.The model formulations have only been tested on a case basis, and although strong signals have been found, this does not imply improved meteorological performance of the model.In particular, testing over longer periods including all seasons was not conducted that time.Furthermore, the interactions between aerosols and the cloud ice phase are not in a state where improvements would be expected.Therefore, it is necessary to mention that it is too early to make conclusions about the improvement of precipitation forecasting by implementation of the indirect aerosol effects, because of large uncertainties in parameterisation of the aerosol-cloud microphysics processes (especially for ice nucleation) and due to adjustments of such effects indirectly in NWP model parameters and constants (retuning of them after implementation of the aerosol feedbacks is needed).More investigations, further improvements and evaluations are needed for aerosol in-direct effects and aerosol-cloud microphysics schemes in the model.Recently such evaluation studies are realised within the CarboNord project for monthly and annual validation studies and will be published separately.

Urban meteorology and environment prediction and assessments
The analysis of urban boundary layer (UBL) for metropolitan areas of the megacity of Paris (population of more than 10 million) and the growing medium-size city of Bilbao (population of 1 million) placed over semi-flat and coastalcomplex terrains, respectively, was performed employing the Enviro-HIRLAM model.In particular, the (1) evaluation of the model performance coupled with the urban module for different types of terrain and sizes of cities and (2) estima- tion of urban heat island (UHI) development over selected urban areas and surroundings were done.The Enviro-HIRLAM simulations were performed for nested domains with horizontal resolutions of 15, 5 and 2.5 km and for selected periods in July 2009.The meteorological boundary conditions were provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) every 3 h.The model was employed in two modes.The first mode is the control (CTRL) run.The second mode is the urban (URB) run -e.g.coupled with the BEP (Martilli et al., 2002) module and anthropogenic heat fluxes (AHFs) from the Large scale Urban Consumption of energY (LUCY) model (Allen et al., 2010).Extracted AHFs were 60 and 40 W m −2 for the Paris and Bilbao metropolitan areas, respectively.For the URB run at the finest resolution, the Paris and Bilbao urban areas were represented by 220 and 16 urban cells, respectively (Fig. 9; adapted from González-Aparicio et al., 2010).In each grid cell, BEP parameterises the flux exchange between the urban surface and the atmosphere depending on combination of different urban districts, e.g.residential, low and high buildings, industrial and commercial.
The statistical analysis showed that the urban simulation had a reduced bias with respect to observations than the control simulations.For Paris, on a monthly basis, the correlations for air temperature were higher for the URB compared to CTRL run, and results improved up to 10 % on a diurnal cycle (with a maximum of 0.83 at 08:00 UTC).The correlations were slightly lower (down to 0.5) at early morning hours and slightly higher (up to 0.8) during afternoon and nighttime.Moreover, correlations at suburban and urban stations were similar to correlations at rural stations (see Fig. 10a).Analysis for Bilbao (González-Aparicio et al., 2013) showed similar performance of the model for both runs: with correlation for air temperature about 0.85 and 0.88 for summer and winter, respectively.For the specific humid-ity, it was 0.75 and 0.92.For the wind speed, the highest value (0.8) is in summer, and during winter it decreased to 0.6 (0.4) near the coast (inland) stations.
The results of simulations for two selected cities showed that the model reproduced well the mesoscale processes at the regional scale for inland winds over Paris and landsea breeze interactions over Bilbao.For selected locations (e.g.coastal versus inland sites), the bias between the observations and simulations was higher over Bilbao (maritime) than over Paris (continental).Although hydrostaticity of the model over a complex terrain is a limitation, but a sensitivity test over Bilbao showed that at 2.5 km optimal resolution it is possible at the same time to obtain satisfactory reproducibility of the large-scale processes and to explore the urban effects at finer scales.
The UHI development was also for short-term periods (here, for Paris -28 July 2009; for Bilbao -15 July 2009) with calm and anticyclonic conditions.For Paris, three different locations were considered: urban (LHVP), suburban (SIRTA) and rural (CHARTRES) stations (Fig. 10a).As seen, the UHI was fully developed at 04:00 UTC with an air temperature anomaly of 2.2 • C (LHVP) and 0.6 • C (SIRTA).It started at midnight and expanded, covering an area of about 2000 km 2 .The heat island was retained until 11:00 UTC, but during the daytime (e.g.11:00-17:00 UTC) the effect disappeared due to the contribution of incoming solar radiation.At CHARTRES, this effect (0.2 • C) was almost negligible.Both the wind speed and relative humidity were also affected by the urban area: at LHVP, the wind speed was reduced by a maximum 3.5 m s −1 at 06:00 UTC, and the relative humidity was down to 15 % under developing UHI.At SIRTA, the change in wind speed was down to 0.7 m s −1 , and at CHARTRES the changes in wind speed and relative humidity were almost negligible.For Bilbao, the model showed that for breezes from northern directions, the impact of urban area on local flow dynamics is inhibited; however, for breezes from southern directions, the urban effect had appeared.For example, on 15 July 2009, the UHI was developed during night-morning hours (e.g.23:00-09:00 UTC) with a maximum up to 1 • C, and the heat island expanded covering area of about 130 km 2 .In addition, González-Aparicio et al. ( 2013) showed that the UHI intensity is lower in winter compared with summer; underlying that dominating factor is the surface heating during daytime, which is higher in summer than in winter.
As medium-size cities are under continuous development, future impacts of urbanisation are expected to become more significant.Several different scenarios of urban development were tested for Bilbao (González-Aparicio et al., 2014).Enviro-HIRLAM model runs showed that under calm conditions during summer and winter, the UHI could reach up to 2.2 • C, covering an area of about 400 km 2 when the city is doubled in size or doubled in AHF.When the city is tripled in size, the UHI could reach up to 3 • C with urban island expansion up to 550 km 2 (Fig. 10c).Analysis of UHI for Bilbao (e.g.triple-size city scenario) versus current UHI over Paris showed similar intensity of up to 3 • C, and UHI boundaries were different; e.g. for Paris, they were was 4 times larger.Such differences can be explained by different cities' sizes, morphologies and characteristic AHFs.

Pollen forecasting
Among air-pollinated allergens, birch pollen is one of the most important for the population group suffering allergic diseases.The number of allergic patients sensitive to birch pollen is assessed as 20 % of European population (WHO, 2003;Linneberg, 2011) and this number is constantly increasing.In particular, in Denmark, the number of allergic patients has increased twice over the past few decades (Lin-neberg, 2011).These facts demonstrate the importance of operational birch pollen forecasting for the European population especially during the spring season.Currently, birch pollen is presented as a biological air pollutant in different NWP and ACT models such as SILAM (Finland), COSMO-ART (Germany, Switzerland, Austria), CHIMERE (France), Enviro-HIRLAM, DEHM (Denmark) and others.The pollen emissions are strongly dependent on meteorology, so it is advantageous to simulate and forecast pollen pollution episodes by online coupled meteorology-air pollution models since all necessary meteorological fields are available at each model time step.
Original developments of the dynamical Enviro-HIRLAM-based operational modelling system for the birch pollen forecasting in Denmark (called Env-POLL) were started in 2006 (Rasmussen et al., 2006;Mahura et al., 2006b) including previously developed statistical methods (Rasmussen, 2002), modelling of elevated concentrations episodes, analysis of spatiotemporal and diurnal cycle variabilities, contribution of remote source regions into pollen levels, improvements in emissions and parameterisations, etc. (Mahura et al., 2007b(Mahura et al., , 2009(Mahura et al., , 2010a)).The most recent developments are shown in Kurganskiy et al. (2015) with a revised general scheme of input and output of the Enviro-HIRLAM birch pollen forecasting system presented in Fig. 11.The input includes the meteorological initial/boundary conditions (IC/BC) obtained from the IFS model system, birch forest fraction map, phenological data, i.e. temperature sum thresholds for the start of flowering (Sofiev et al., 2013), and the accumulated total number of birch pollen particles emitted from a unit area during the pollinating season.
The forecasting of birch pollen concentrations requires information/data on the spatial birch tree distributions, characteristics of pollen release, its atmospheric transport and dispersion, its deposition due to gravitational settling and  wet deposition, i.e. scavenging by precipitation.Birch pollen emissions are fully dependent on temporal and spatial variability of meteorological conditions.The emission module (Sofiev et al., 2013) includes the following parameters affecting the pollen release: 2 m air temperature and relative humidity, 10 m wind speed and accumulated precipitation.The atmospheric transport is handled in the same way as for aerosols (see Sect. 2.8).Dry deposition of birch pollen particles in the atmosphere is represented by gravitational settling (Seinfeld and Pandis, 2006), whereas dry deposition due to interactions of particles with the surface can be neglected according to Sofiev et al. (2006).The wet deposition scheme distinguishes between in-cloud (Stier et al., 2005) and belowcloud scavenging (Baklanov and Sørensen, 2001).The out-put in terms of birch pollen forecasting, for analysis, contains 2-D fields of the birch pollen concentration at the lowest vertical model level.The modelling domain has 15 km horizontal resolution with 154 and 148 grid points along longitude and latitude, correspondingly.The domain covers the main European part and is centred around Denmark.
As examples for the birch pollen season of 2006, the model results were compared with observations for two Danish sites: Copenhagen and Viborg (Fig. 13).This year was dominated by a relatively cold spring over large areas of Europe followed by rapid warming and little/no rain.It caused a short but intensive birch pollen season with long-range transport episodes before the local flowering start and thereby emissions.The evaluation for both modelled and observed birch pollen concentrations showed extremely high values (daily averages about and even more than 1000 grains m −3 ) during the 5-10 May 2006 episode for Copenhagen and 5-8 May 2006 episode for Viborg.The extremely high birch pollen concentrations over Denmark are also visible in Fig. 12b.
According to Sofiev et al. (2011) and Siljamo et al. (2013) the following criteria can be used for assessment of birch pollen concentration forecasting: model accuracy (MA), hit rate (HR), false alarm ratio (FAR), probability of false detection (POFD) and odds ratio (OR).All of the criteria are calculated using four parameters obtained by assessment of the number of low and high modelled versus observed birch pollen concentrations (C) relative to a threshold value N th = 50 grains m −3 (i.e.C ≥ N th for high-concentration and C < N th for low-concentration days).The threshold has been chosen since most of the pollen-allergy-sensitive population might start suffering from allergic reactions when daily mean birch pollen concentration C is greater than or equal to N th in the air (Jantunen et al., 2012).
The results of statistical analysis showed high MA for both Danish stations (0.95 for Copenhagen and 0.84 for Viborg, 0.9 on average).Prediction of elevated/top concentrations (HR values) by the model was assessed as 0.93 for Copenhagen and 0.58 for Viborg.The FAR values indicated that the probability to get an incorrect top model concentration was 0.07 and 0.42 for Copenhagen and Viborg, respectively.The POFD criterion showed low probability to get high modelled concentrations for observed low-concentration days (0.02 for Copenhagen and 0.18 for Viborg).Finally, the OR indicated that the likelihood of getting "high" day concentration instead of a "low" (if the model prediction is "high") was 42 and 3.26 times higher for Copenhagen and Viborg, respectively.In other words, the OR values show the ratio between HR and POFD.As it is seen from the OR values provided above, a fraction of the correct forecasts is prevailing for both Danish stations in this study.
It was found that, compared to observations, the modelled results reflected the general shape of changes in pollen concentration during the episode studied for both Danish stations: Copenhagen and Viborg.As it is also seen in Fig. 13 that the model reproduces the magnitude of birch pollen concentrations for the peak period of the season in comparison with observations.However, some overestimation of the modelled concentration is visible for both stations at the end of the season.It can be explained by contribution due to longrange atmospheric transport of pollen from other remote regions, presumably from those located more northerly than Denmark and where the pollen season starts and ends later relative to the Danish sites.

Chemical weather forecasting and air pollution applications
Validation and sensitivity tests (on examples of case studies and short time episodes) of the online versus offline integrated versions of Enviro-HIRLAM (Korsholm et al., 2008) showed that the online coupling improved the results.Different parts of the model were evaluated versus the ETEX-1 experiment, the Chernobyl accident and Paris MEGAPOLI campaigns (summer 2009) datasets and showed that the model had performed reasonably well (Korsholm, 2009;Korsholm et al., 2009Korsholm et al., , 2010;;Sokhi et al., 2017).Online versus offline coupled simulations for the ETEX-1 release showed that the offline coupling interval increase leads to considerable error and a false peak (not found in the observations), which almost disappears in the online version that resolves mesoscale influences during atmospheric transport and plume development (Korsholm et al., 2009).Further studies (Korsholm et al., 2010) of urban aerosol effects on the atmospheric composition showed that aerosol feedbacks through the first and second indirect effects induce large changes in chemical composition, in particular nitrogen dioxide, in the case of convective clouds and little precipitation.For the Paris campaign, on diurnal cycle variability, the ozone concentration patterns showed dependencies on meteorological parameters and especially at urban scale runs (Mahura et al., 2010b).
To perform further analysis of online coupling and feedback effects on atmospheric pollution forecasting, the year 2010 was selected (for details, see Sect.2.5).Nuterman et al. (2013)    and Spain (see Fig. 14a).The model runs were performed for the entire month of July 2010 with a 7-day spin-up in June.
Figure 14b shows correlation coefficients on a diurnal cycle for PM 2.5 concentrations at selected measurement sites.In general, it shows fairly good positive correlations (more than +0.3), except for several Spanish stations (such as ES1938A at daytime and ES1974A at nighttime).On the monthly-based evaluation, the model predicts well PM 2.5 day-to-day variability but always has negative bias (Fig. 15).This underprediction is due to several reasons: (i) aerosol microphysics without secondary organic aerosols; (ii) lack of partitioning of ammonium nitrate; (iii) rough model resolution, which still cannot capture small-scale effects like complex orography and urbanised regions (in particular, due to lack of fine-resolution emissions from anthropogenic sources, like urban traffic).For instance, the model shows a negative bias of PM 2.5 during daytime at the Danish urban station (Fig. 15a).It is apparently due to rough model resolution in the considered runs.It was also found that PM 2.5 values are very influenced by changes in atmospheric stability conditions, which are difficult to predict accurately in many NWP models.This can be observed from correlation coefficient decrease at stations during nighttime (at 03:00 UTC) or from underestimation of elevated concentrations.In spite of these issues, the model can well reproduce the diurnal cycle of aerosols at different sites, e.g.urban (Fig. 15a), coastal and rural (Fig. 15b), and shows good overall performance.
Further on-going developments of the Enviro-HIRLAM modelling system for atmospheric composition applications are realised within the FP7 MarcoPolo and NordForsk Car-boNord projects.The Enviro-HIRLAM downscaling from regional-to urban-scale modelling is realised in MarcoPolo for the east China region and the largest metropolitan agglomerations in China (Mahura et al., 2016) with a focus on providing services on meteorology and atmospheric composition (with a focus on aerosols).The Northern Hemisphere low-resolution modelling in a long-term mode is realised in  CarboNord with a focus on evaluation of black carbon as well as higher-resolution modelling over the European domain in a short-term mode with a focus on feedback mechanism evaluation (Nuterman et al., 2015;Kurganskiy et al., 2016).

Further discussions
Several types of the above-described and previously published applications of the Enviro-HIRLAM for meteorological, environmental and climate forecasting, and assessment studies were tested and demonstrated.Different applications of Enviro-HIRLAM (with downscaling from hemispheric -regional -subregional -urban scales) were realised for different geographical regions and countries including European countries such as Denmark, Lithuania, France, Spain, Ukraine, Russia, the Netherlands and Turkey, as well as for others (Kazakhstan, China and Arctic regions).
It is clear that the seamless/online integrated modelling approach realised in Enviro-HIRLAM is a perspective and state-of-the-art way for future single-atmosphere modelling systems, providing advantages for all three communities: meteorological modelling including NWP, air quality (AQ) modelling including CWF, and climate modelling.
However, there is not necessarily one configuration of the integrated online modelling approach/ system suitable for all www.geosci-model-dev.net/10/2971/2017/Geosci.Model Dev., 10, 2971-2999, 2017 communities, and that should be further investigated with practical needs for area applications, approaches to coupling and computing resource usage.In particular, based on previous studies and above-shown examples, the following could be recommended for the considered applications: -For AQ, online coupling improves air quality forecasts, especially with full chemistry and aerosol feedbacks effects included.
-For NWP, gas chemistry is not critical and can be simplified (or omitted), but aerosol feedbacks are important for radiation and precipitation and especially for heavy polluted episodes and in urban areas.
-For pollen forecast, online coupling improves pollen emission parameterisation and correspondingly modelling of concentration and deposition; however, the feedbacks are not so important.The chemistry is not considered yet, but interaction with allergens would be interesting to study in future (this has not been done yet).
-For climate studies, it is suitable only for understanding the feedback mechanisms but too expensive computationally for climate timescale runs (the model had been used usually for 1-year runs); chemistry is important, as there is a need for it to be optimised and simplified.
It should also be mentioned that the considered evaluations of Enviro-HIRLAM were done only for some elements (e.g. the coupling interval) in the previous analysis and main conclusions about the improvements were provided just for these.Other feedback mechanisms, especially for aerosol-cloud interactions, were analysed mostly as sensitivity studies or evaluated for short-term episodes.In particular, the STRACO cloud scheme contains fairly simplified cloud microphysics (heavily parameterised).Hence, tuning is essential for the overall performance of the model when it comes to precipitation and cloud physical properties.

Conclusions
In this paper, we have provided a comprehensive description of the Environment -High Resolution Limited Area Model (Enviro-HIRLAM), which is developed as a fully online coupled/integrated numerical weather prediction and atmospheric chemical transport modelling system for research and forecasting of joint meteorological, chemical and biological weather.
Possible applications of the modelling system can include chemical weather forecasting, air quality and chemical composition for short-and long-term impact assessments on population and environment, multiscale weather forecasting (e.g. on regional and subregional scales, in urban areas and severe weather events), pollen and road weather condition fore-casting, climate change forcing modelling, studies of climate change effects on atmospheric pollution on different scales, weather modification and geoengineering methods, forest fires and volcano eruptions, dust storms, nuclear explosion consequences and other emergency preparedness modelling.
Comprehensive online modelling systems, like Enviro-HIRLAM, built originally for research purposes and including all important mechanisms of interactions, will help to understand the importance of different physical-chemical processes and interactions as well as to create specific model configurations that are tailored for their respective purposes.
Multiple-episode studies with the Enviro-HIRLAM model demonstrated the importance of including the meteorology and chemistry (especially aerosol) interactions in online coupled models.However, there is no unique integrated online modelling system configuration which is best suited for all communities.
Highlighting a number of previous investigations, we show that Enviro-HIRLAM has already been used for a host of different applications ranging from pollen forecasting to numerical weather prediction.
It should be stressed that there are still main gaps remaining in the understanding of several processes such as (i) aerosol-cloud interactions (still poorly represented); (ii) data assimilation in online models (still to be developed to avoid overspecification and opposite cancelling effects); and (iii) model evaluation for online models (which needs more process data and long-term measurements, and a test bed).

Figure 1 .
Figure 1.General scheme of international collaboration, research and development, technical support and science education for the online integrated Enviro-HIRLAM: "Environment -High Resolution Limited Area Model".

Figure 3 .
Figure 3. (a) The typical SW spectrum used for calculating average SW aerosol optical properties.(b) The spectral weights used for calculating average LW aerosol optical properties.

Figure 4 .
Figure 4. General scheme of the building effect parameterisation (BEP) module for the Enviro-HIRLAM model urbanisation with a structure of the sub-routine conception (adapted from Mahura et al., 2010b).

Figure 5 .
Figure 5. Example of the simulated PM 2.5 concentration over Europe on 17 July 2010 with horizontal resolution of 16 km.

Figure 9 .
Figure 9. Urban district classification based on urban zoning data for the (a) Bilbao and (b) Paris metropolitan areas, including the residential area (ReD), low and high building districts (LBDs and HBDs, respectively) and industrial and commercial districts (ICDs).Spatial distribution of urban districts (HBDs -high buildings, RDs -residential, ICDs -industrial commercial and CCs -city centre districts) for the Paris metropolitan area within the P01 modelling domain (partly adopted from González-Aparicio et al., 2014) is shown.

Figure 10 .
Figure 10.Difference plots for the air temperature at 2 m between outputs of the URB (urbanised BEP plus AHF) and CTRL (non-urbanised) Enviro-HIRLAM model under calm conditions during summer 2009 for the (a) Paris metropolitan area and for the Bilbao metropolitan area (b) in its current size of the city and (c) under a scenario tripling the size of the city.

Figure 12 .
Figure 12.(a) Birch forest fraction map; (b) example of the simulated birch pollen concentration in the modelling domain on 6 May 2006 at 18:00 UTC.

Figure 15 .
Figure 15.Error-bar concentrations (µg m −3 ?) on the diurnal cycle for AirBase observations versus Enviro-HIRLAM modelling results: (a) Danish urban station and (b) German rural station.The top right corner indicates maximum and average correlation coefficients for the station as well as total number of analysed observation samples.Green numbers along the x axis indicate number of observation samples per time slice.