The On-Line Integrated Mesoscale Chemistry Model BOLCHEM

: This work presents the on-line coupled meteorology–chemistry transport model BOLCHEM, based on the hydrostatic meteorological BOLAM model, the gas chemistry module SAPRC90, and the aerosol dynamic module AERO3. It includes parameterizations to describe natural source emissions, dry and wet removal processes, as well as the transport and dispersion of air pollutants. The equations for different processes are solved on the same grid during the same integration step, by means of a time-split scheme. This paper describes the model and its performance at horizontal resolution of 0.2 ◦ × 0.2 ◦ over Europe and 0.1 ◦ × 0.1 ◦ in a nested conﬁguration over Italy, for one year run (December 2009–November 2010). The model has been evaluated against the AIRBASE data of the European Environmental Agency. The basic statistics for higher resolution simulations of O 3 , NO 2 and particulate matter concentrations (PM 2.5 and PM 10 ) have been compared with those from Copernicus Atmosphere Monitoring Service (CAMS) ensemble median. In summer, for O 3 we found a correlation coefﬁcient R of 0.72 and mean bias of 2.15 over European domain and a correlation coefﬁcient R of 0.67 and mean bias of 2.36 over Italian domain. PM 10 and PM 2.5 are better reproduced in the winter, the latter with a correlation coefﬁcient R of 0.66 and the mean bias MB of 0.35 over Italian domain.


Introduction
For many years, Meteorology and Air Quality (AQ) have been developed as separate disciplines. From the modeling point of view, this led to the development of separate Meteorological and Chemical Transport Models (CTM), the latter using as input data the output of the former. However, over the last decades, it was recognized that the interaction among processes in the atmosphere could not be confined within those two domains. Actually, some of the missing couplings in the offline representation of the atmosphere proved to be fairly important, not only for climate modeling where the impact of greenhouse gases and aerosol on radiation and clouds is of paramount importance in the determination of the global (and local) budgets, but also for short time forecasts of both meteorological and composition quantities, in which aerosols play a major role (see [1] for a review).
Despite this, the development of integrated models of atmospheric dynamics and composition is relatively recent: for instance, the first stable version of WRF-Chem, one of the most widespread models, has been released only in 2005 [2]. In Baklanov et al. [1] the status and a short model description of the more used online coupled models in Europe is presented. Among others, we recall a.e., the COSMO-ART model [3], Enviro-HIRLAM model [4], the Me-tUM [5,6]. The online coupled models involve a strong coupling of the two "realms" with raising scientific and computational complexities. For these reasons, CTMs, i.e., offline models that typically use the output of a meteorological model to compute transport, diffusion, and transformation of pollutants, are still widely used (see, e.g., Copernicus Atmospheric Monitoring Service-Europe, https://atmosphere.copernicus.eu).
In recent years, large efforts have been undertaken to develop not only integrated models, but also strategies and frameworks for online integrated modeling (see, e.g., Baklanov et al. [7]) and their evaluation, as, for instance, the joint European-North American Air Quality Model Evaluation International Initiative [8,9].
BOLCHEM is an online integrated mesoscale meteorology-chemistry model where, as defined in Baklanov et al. [1], equations for different atmospheric processes are solved "simultaneously" (by means of a time split scheme) over the same grid using one main time step for integration. Although in the online integrated models it is possible to include interactions between different components (typically "feedback"), at the present stage of development, in the BOLCHEM model, the chemistry and meteorological components are one-way coupled (defined as "online-access" according to Baklanov et al. [1]). The processes that constitute the meteorological component (flow dynamics, thermodynamics, radiation, surface and other physical processes) are based on the BOLAM model, which is a meteorological model developed in Italy, at ISAC-CNR. BOLCHEM development started in 2002, making it one of the first projects for the development of online integrated models in Europe. The first version of BOLCHEM included only the gas-phase chemical processes. BOLCHEM was inserted in the NetFAM/COST 728 database among the "few" (at that time) online models (see Baklanov et al. [10]). From 2005 to 2009, BOLCHEM participated in the GEMS Project [11], the Regional Air Quality (RAQ) subproject with a version that included gas chemistry and related processes [12,13], and to MACC-Modelling Atmospheric Composition and Climate. A first study with BOLCHEM can be found in Mircea et al. [14]. In the CityZen project, BOLCHEM contributed to the activities on regional climate [15][16][17], including aerosol and its interaction with radiation. The model, with the inclusion of both gas and aerosol chemistry, was applied for specific air quality studies, mainly in the frame of INTERREG ERESIA, INTERREG CESAPO, and MED POSEIDON projects [18][19][20].
This paper presents the first full description of the meteorological and chemical modules of the model (Section 2) and its performances on the European and Italian domain (Section 3).

Model Description
Every existing integrated meteorological and chemical transport modeling system (Met-Chem Model) is built on top of a pre-existing meteorological model. This is reflected in the model description. Here we start the description keeping in mind the online-integration framework, outlining the model description based on processes: turbulence (including dispersion of tracers), interface (including dry deposition), radiation (including primary effect), cloud (including wet removal), and so on. BOLCHEM is based on the hydrostatic meteorological model BOLAM [21], the gas-phase chemical mechanism SAPRC90 [22] and the aerosol model AERO3 [23,24] and includes a number of parameterizations to describe transport and dispersion, natural emission sources, dry and wet removal processes and their interactions. These different components have been intimately connected by unifying as much as possible the process involved. All the meteorological and chemical fields (gas concentration, number and mass concentration of aerosol species, virtual temperature, specific humidity, . . . ) are available at every time step of integration (see Section 2.7). This allows a consistent treatment of the different processes, as advection, turbulence, radiation, for both meteorological and chemical quantities.
This section describes in detail the model components, the modeling system configuration and applications.

Boundary Layer and Vertical Diffusion
The surface layer (SL) is modeled according to the classical Monin-Obukhov similarity theory [25]. The Businger (see [26]) stability functions are used in the unstable SL, while Holtslag [27] functions apply to the stable case. The roughness length over land, initially defined depending on vegetation and sub-grid orographic variance, is modified as a function also of snow coverage conditions. Over the sea, a Charnock roughness [28] representation is introduced for computing momentum fluxes. It takes into account the dependence of wave height on the surface wind speed, while roughness lengths for temperature and humidity in stable and unstable conditions are defined according to [29]. The mixing length (ML) based turbulence closure, widely used to compute the PBL fluxes for atmospheric modeling (see, for instance, [30]) is applied to model the turbulent vertical diffusion of momentum, potential temperature and specific humidity in the free atmosphere. The turbulence closure is of order 1.5, in which the turbulent kinetic energy (TKE) equation is integrated over time [31]. To take into account buoyancy effects in case of saturated atmosphere, the ML depends on the Richardson number based on equivalent potential temperature. In the unstable case, a modified version of the Bougeault and Lacarrere [32] ML is applied while, in the stable case, a modified Blackadar [33] formulation is used. Finally, the dissipated TKE is fed back into resolved temperature (frictional heating).
Diffusion is split into vertical and horizontal. The vertical one is performed using the vertical diffusion coefficients for all of the atmospheric components (temperature, humidity, gas and aerosol species). This avoids the need to compute extra diagnostics variables typically used in CTM (e.g., the mixing height). Furthermore, the calculation of all the derived quantities (i.e., U 10 ) is straightforward in a coupled model where all the quantities are available online/runtime.

Convection and Precipitation
The sub-grid convective precipitation is treated in BOLCHEM following the Kain-Fritsch (KF) parameterization scheme [34,35] that has shown considerable success in simulating the development and evolution of convection under a variety of atmospheric environments [36,37]. The parameterization has been completely re-coded and some modification has been introduced in the original code. Liquid water static energy (instead of the Bolton approximation of equivalent potential temperature) is used as a thermodynamic conserved quantity. Moreover, additional modifications have been introduced with respect to the Kain [35] version, regarding the dependency of downdraft on ambient relative humidity, and the precipitation rate. The cloud depth threshold establishing the onset of shallow convection has also been increased. The above changes tend to slightly reduce, on average, the temperature at low tropospheric levels around and below cloud base, hence stabilizing a little more efficiently the lower troposphere. This has also the effect of reducing to some extent the intensity of small-scale cyclogenesis in the presence of convection.
The scheme describing the physical processes of stratiform precipitation is based on a single-moment, bulk microphysics scheme with two ice precipitation categories (snow and graupel). The spectral properties of clouds and hydrometeors are described assuming generalized gamma function distributions. The main processes described by the microphysical scheme are: nucleation of cloud water (cw) and of cloud ice (ci), condensation and evaporation of cw, freezing of cw, sublimation and melting of ci, auto-conversion of cw and of ci into rain and precipitating ice categories, sublimation of snow and graupel, evaporation of rain, melting/freezing of hydrometeors, hydrometeor and cloud interactions (collection/accretion/riming), computation of terminal fall speeds and fall process, using a conservative backward-upstream integration scheme, thermodynamic feedback based on enthalpy conservation. Aerosols and gases removal by clouds and precipitation is treated as in the EMEP/MSC-W model [38]. The online availability of the meteorological fields allows a better computation of the effects of spatial and temporal variability of clouds and precipitation on pollutant removal, since the evolution of single-cell storms have a typical lifespan of less than one hour that cannot be properly captured by offline models.

Gas Chemistry
The atmospheric gas-phase chemistry is simulated with the SAPRC90 photochemical mechanism, as described in [22]. It is implemented using the chemistry preprocessor Flexible Chemical Mechanism (FCM) [39]. This mechanism considers 35 chemical species that undergo 131 reactions. It is designed for an efficient representation of complex ambient mixtures with the purpose of assessing differences in atmospheric impacts of individual VOCs. Although it is recognized that chemical processes involving VOCs and NOx lead to the generation of ozone and other secondary pollutants in the presence of sunlight [40,41], the role of VOCs in ozone and secondary organic aerosol formation and evolution is a subject in continuous progress. Different mechanisms use different model species to represent how organic compounds react [42]. SAPRC90 mechanism includes reactions for a large number of VOCs that have been evaluated against environmental chamber data. The SAPRC90 photolysis rates used in photolytic reactions are calculated by the model as a function of the simulated radiation (see Section 2.5). The mechanism includes the gas-phase reactions of SO 2 , but not heterogeneous reactions. SAPRC90 was extended to describe the formation of condensable organic products and the oxidation of SO 2 in aqueous, following the scheme reported in [43].

Aerosol Dynamics
Chemical composition of aerosol is simulated with the aerosol dynamic model AERO3 [23,24] coupled with the inorganic, thermodynamic equilibrium model ISORROPIA [44], and with the secondary organic aerosol model SORGAM [45]. The aerosol dynamic model is coupled with the gas-phase chemical model allowing the use of updated gas-phase precursor concentrations. The implemented aerosol dynamic module follows the so-called modal approach proposed by Whitby [46], in which the population of particles is described by a superposition of log-normal distributions called modes. The aerosol population is described by three modes (Aitken, accumulation, coarse), the aerosol particles are formed by nucleation and growth by condensation and coagulation, both intra-and inter-modal. Gas/particle mass transfer continuously modifies the chemical composition of aerosol. The considered aerosols are: primary unspeciated anthropogenic aerosol and water, primary organic aerosol, primary elemental carbon, soil and sea salt, secondary biogenic and anthropogenic organic aerosol, ammonium, nitrate, and sulfate. The amount of ammonium, nitrate, sulfate and water contained in the aerosol is calculated with ISORROPIA [44] and depends on available concentrations of ammonia, nitric acid, SO 2 and on relative humidity and temperature. The amount of secondary organic aerosol (SOA), both biogenic and anthropogenic, is calculated using the absorptive partitioning model of Pankow [47], extended by Odum et al. [48], and the concentrations of reacted ROG's, the total (gas + aerosol) concentration of each semi-volatile organic compound and the concentrations of primary and secondary organic aerosols.

Radiation
The atmospheric radiation is computed with a combined application of the Ritter and Geleyn [49] scheme (hereafter referred to as RG) and the Morcrette [50], Morcrette et al. [51] scheme, (referred to as ECMWF), the latter including the Rapid Radiative Transfer Model (RRTM) [52] and the Tegen et al. [53] aerosol climatology. The RG scheme, with the option of maximum cloud coverage, is called approximately every 0.5 h of simulated time and has been modified to take into account the explicit cloud water and cloud ice concentrations predicted by BOLCHEM (in the original RG scheme, clouds were parameterized from specific humidity). The ECMWF scheme, being much more computationally expensive than the RG scheme, is computed only every 1.5 h of simulated time at alternate horizontal grid points, in order to save computational time. The difference between the surface and internal radiative fluxes of the ECMFW and RG schemes is then interpolated in time and space, and it is used to correct the fluxes computed by the RG scheme. Concerning photolytic reactions, the computation of photolysis rates requires the use of algorithms that make use of spectral radiative and molecular properties (absorption cross-sections and photo-dissociation quantum yields) [54] that are rather expensive to be computed run-time.
One of the typical alternative approach (see [1], Table 8), i.e., is to derive the photolytic rates from a lookup table that contains values previously computed in clear-sky conditions as a function of the solar zenith angle. This approach is used in BOLCHEM, with photolytic rates produced by running the FCM over 11 h angles. The clear-sky values are then modified locally (in every grid cell) according to the ratio between the actual short-wave net flux and the same quantity in clear sky conditions as computed by the model, both available runtime. This allows for a 3D spatial distribution of photolysis rates that accounts for the actual (modeled) local radiation fluxes. The correction factor can exceed 1.0 above clouds and is reduced within and below clouds.
Because the impact of clouds on the actinic flux is nearly independent on the wavelength [55] the correction is applied without reference to specific spectral properties of involved molecules, and it is assumed to capture the overall modification of the photolysis activity without adding heavy computation. A similar approach was adopted in other models, e.g., COSMO-ART [3].

Surface Fluxes
Surface fluxes (sources and sinks) are the most critical input for a Met-Chem. Moisture and heat influence atmospheric dynamics while anthropogenic and natural emissions are the source of gas and particulate pollutants. On the other hand, removal by dry deposition can be an important sink for them. In the following, submodels of surface fluxes for the different components are described. While anthropogenic contribution can be only included from available databases, natural component sources (heat, moisture, gases and aerosols) couples different components such as radiation, turbulence, soil type and moisture content. Different parameterizations are used for different processes. According to the online approach, all of the modules used for the calculation of surface fluxes of the various components share all the available information.

Anthropogenic Emissions
The anthropogenic emission data are based on European Database or National Database. In the present work, we use the TNO-MACC-III emission database for the year 2010, that is an update from the TNO-MACC-II dataset [56]. It provides gridded emissions on the European domain, at 0.125 • × 0.0625 • longitude-latitude resolution, for each year, country, sector and source type. Additionally, a classification in area and point sources is provided without any information on emission height and other parameters, for individual point sources. Eight key families of pollutant are considered: methane (CH 4 ), carbon monoxide (CO), ammonia (NH 3 ), non-methane volatile compounds (NMVOC), nitrogen oxides (NO x ), particulate matter (PM 2.5 and PM 10 ), sulfur dioxide (SO 2 ). The annual emissions are then split in time to hourly resolution and NMVOC and PM are speciated on the basis of emission sectors (SNAP) and seasonal activity [57], according to the chemical module SAPRC90 and the aerosol module AERO3. Since within one SNAP, different types of sources occur, the emissions can have different height. We have used the following recipe to allocate emissions on vertical levels: area sources at the lowest layer of the model, point source vertically allocated on one of the three levels (0 m, 50 m, 150 m) depending on the typical stack emission for different SNAP.

Soil and Vegetation Model: Moisture and Heat
The meteorological component BOLAM includes an original soil model with fixed lower boundary conditions, where heat and water vertical transfer are computed at the interfaces of several soil layers with depths ranging typically from a few cm to more than 1 m, increasing downward. Vegetation effects at the surface (transpiration and interception of precipitation) and in the soil (extraction of water by roots depending on wilting conditions), taking into account different soil types and physical parameters, are considered. The soil model also includes a treatment of freezing and melting processes of the water content. At the surface, the evolution of the snow cover is computed considering snow accumulation and melting, with a single-layer snow mantle model. A surface skin temperature is defined by imposing zero-net flux divergence of heat at the soil-atmosphere interface. Albedo and emissivity variations are also computed as a function of the uppermost soil water content. The sea surface temperature is predicted using a slab ocean model, where latent and sensible heat fluxes, and radiation contributions, are taken into account. The soil and vegetation dataset is obtained from the Global Land Cover Facility of the University of Maryland at 1/120 degree resolution [58], a classification from AVHRR satellite and free available.

Biogenic Gas Emissions
The surface flux of visible radiation and skin temperature, calculated run time by the BOLCHEM model at every time step, are used to estimate the biogenic gas emission fluxes. Biogenic emissions are grouped into three categories: isoprene (C 5 H 8 ), monoterpenes (C 10 H 16 ) and OVOCs (C x H y O z ). Following Guenther et al. [59], Simpson et al. [60], the biogenic emission flux is given by the product of the emission potential, the foliar biomass density and an environmental correction factor representing the effects of temperature, and in some case solar radiation, on emissions. The emission potential and the foliar biomass density are ecosystem-dependent. BOLCHEM uses an independent database, having 1 km spatial resolution [61]: a non-canopy approach is adopted and the emission potential is calculated at branch-levels. For monoterpene, two separated emission potentials are provided. Depending on the plant specie, the monoterpene emission flux is driven by temperature only, or by light and temperature. The correction factor is calculated following Guenther et al. [62]. For isoprene, as for monoterpene light and temperature depending, the factor is the product of two coefficients which account respectively for the effect of leaf temperature and Photosynthetically Active Radiation (PAR). PAR is calculated in BOLCHEM as a function of the surface flux of visible radiation. For monoterpene and OVOCs emission from most plants, the correction factor is a factor depending on temperature only [59,62]. BOLCHEM uses the skin temperature to calculate the correction factor.

Sea Salt
Sea salt emission flux is calculated following Zhang et al. [63]. In this approach, the solute weight fraction of seawater, composed in the model by Cl − , Na + and SO 2− 4 , is represented as a function of RH [64] over the 0.45-0.99 RH range. Then, the density of a sea-salt particle is expressed as a function of the solute weight fraction of seawater. The surface fluxes of sea-salt are first calculated at the reference relative humidity of 80%, using open-ocean source functions that empirically relate the particle number size distribution and the wind speed at 10 m a.s.l. (U 10 ) [65,66]. This function is applicable for U 10 < 20 ms −1 and 0.8 m < r 80 < 10 m, being r 80 the wet particle radius at RH = 80%. The size distribution of sea salt particle flux is then adjusted according to the local relative humidity, using a correction factor, expressed as a function of ambient relative humidity RH. The relative humidity RH and U 10 are calculated by BOLCHEM at every time-step, the latter according to the stability condition (see Section 2.1).

Forest Fire Emissions
Forest fire emissions can be estimated for particular case studies and then included in the BOLCHEM model [18]. The pre-processor prebolchem_fire is described in Pizzigalli [67]. Wildfire emission fluxes are estimated starting from latitude and longitude of area burned obtained from satellite products, as MODIS burned area product, following the method-ology proposed by Seiler and Crutzen [68]. The emission E(x) of specie x is given by the product of the burnt area (A), the fuel load (B), the burning efficiency CE and the emission factor e(x) of species x. B and e(x) are functions of the land cover classification and CE is a function of the tree cover. For each class of vegetation, values of fuel load and burning efficiency are estimated for each class of vegetation, referring to the UMD Global Land Cover Classification (GLCC) [58]. The WRAP approach (WRAP, 2005) is used to modulate the emission fluxes and to estimate the fire emission height. Basically, the fires are classified into size classes based on virtual acreage, then top and bottom emission height (Htop max and Hbot max ) and buoyant efficiency BE size are assigned to every plume class. The total fire emissions are split into the smoldering and flaming part. The smoldering fraction is calculated as a function of the BE size and it is daily modulated as a function of the BE size and buoyant efficiency. The remaining part is distributed between the bottom and the top of the plume following a diurnal profile, with peak emissions during the early afternoon and low emissions during the night.

Dry Deposition
Parameters from the BOLAM surface layer, such as roughness length and stability functions (see Section 2.1) temperature, pressure, soil and vegetation types are directly provided at every time step to the gas dry deposition module. The scheme is based on a resistance analogy approach [69] in which the dry deposition velocities are written as the inverse of the sum of three resistances, Vd = (R a + R b + R c ) −1 , where the three quantities R a , R b and R c are the aerodynamic, quasi-laminar and canopy resistances, respectively. The formulation of R c follows the methodology implemented in Anav et al. [70]. It uses the parameterization of Emberson et al. [71] to account for stomatal conductance, which is calculated for different land cover types as a function of the vegetation type, phenology, and surface parameters as temperature and humidity. The algorithm is adapted to the BOLCHEM soil and vegetation classes, using 16 vegetation and 9 soil types.

Numerical Core
The model prognostic variables are distributed in the vertical on a non-regular Lorenz grid [72], with higher resolution in the atmospheric boundary layer near the lower surface. The vertical coordinate is based on a hybrid, pressure-based coordinate system, in which the pure terrain-following coordinate at the model surface, gradually tends to a pure pressure coordinate with increasing height above the ground. The horizontal discretization is based on a staggered Arakawa-C grid [73], in geographical coordinates (latitude-longitude). The equator can be taken to be any great circle on the Earth in order to minimize grid anisotropy over the limited area of the simulation. The time scheme of BOLAM is splitexplicit, forward-backward for the gravity modes. The main model time step dt is limited by the Courant number constraint for the horizontal advection in the meteorological component. This time step is used for almost all of the processes except in radiation and convection modules, where the computation of the time tendencies of these processes is performed every 30 min and applied to each time step dt, and in the gas chemistry module, where the solver for chemistry uses an algorithm to adapt the time step for which dt represents the maximum. The SAPRC90 equations are solved by the lsode solver with adaptive time-step in which a minimum value is set as an internal parameter. The lateral boundary conditions are applied on a number of grid boundary frames, using a relaxation scheme [74] that efficiently absorbs wave energy, helping in reducing spurious reflection by the lateral boundaries.
The advection scheme presently implemented is the Weighted Average Flux scheme (WAF) [75], which is a second-order implementation of the Godunov method [76] particularly suited to integrate in time the conservation of a scalar quantity. This scheme is a "total variation diminishing" one, and therefore prevents the occurrence of spurious oscillations (see also [77]). The horizontal velocity components are interpolated over the T-points of the Arakawa C grid, advected as any other variable defined on those points, and interpolated back to velocity points (to avoid smoothing, one-dimensional fourth-order interpolation is used).
Conservation of mass of advected and diffused species is easily implemented in an online model. In order to conserve mass exactly in the transport of the scalar parameters involved in the chemical core, a flux-form of the BOLAM advection scheme has been used in which the advected quantities are expressed as volumetric concentrations.
Lateral boundary conditions for meteorology are imposed using a relaxation scheme [74,78], while a sponge layer is adopted at the top of the atmosphere, which is located at about 10 h Pa (middle stratosphere). Boundary concentrations for chemical species are injected into the model by means of an inflow/outflow condition applied to the external grid frame. In addition, to save CPU time, the top level for chemistry can be optionally chosen to be lower than that for meteorology.

Model Configuration and Previous Applications
BOLCHEM has been used in different applications, from experimental forecast over Europe (GEMS) and Italy (MACC) to long term (climatological) studies (CityZen). In these applications, the horizontal resolution ranges from 0.5 • to 0.1 • with a number of vertical levels ranging from 30 to 60. The meteorological component of BOLCHEM uses the hydrostatic approximation and thus the horizontal resolution is limited to about 5-6 km. In order to allow for timesteps as large as possible, horizontal discretization is performed on a rotated, latitude-longitude grid that minimizes the convergence of meridians in the limited area domain. The top of the vertical chemical domain can be selected to be lower than that for meteorology, which is at about 10 h Pa, depending on the case studied. For long range transport simulations and/or for comparison with satellite data, the top can be set to the 200 h Pa or even at the top of meteorological one [79]. For studies at mesoscale the top can be reduced to half of the maximum value of the vertical coordinate. The model timestep is determined by the empiric rule dt 1000 × dlon (s), dlon being the grid distance in degrees. Initial and boundary conditions for meteorology can be taken either from the European Centre for Medium-range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) or from the National Centers for Environmental Prediction (NCEP) Global Forecasting System (GFS), the latter being used typically in forecast mode for real time applications.
As far as the atmospheric composition is concerned, boundary conditions can vary depending on the specific application. In the GEMS and MACC projects, the output from global and regional forecasting systems, respectively, were used. Chemical boundary conditions can be derived by CAMS products (both on global and regional scales). Another possibility is to perform simulations by nesting BOLCHEM into itself. The parent simulation, performed using climatological BC, then provide the boundary conditions for nested simulation. The latter is usually performed at the finer horizontal resolution, but with the same vertical levels as the parent domain. Results of the model used in nesting configuration, with child domain at 0.06 • resolution, are shown in Buccolieri et al. [19] and in Cesari et al. [20]. Anthropogenic emissions must be provided by the user for the specific case considered. In the actual implementation, the used biogenic emissions are based on an inventory generated by NKUA (National and Kapodistrian University of Athens) in the frame of the GEMS project [61]. However, following the same methodology, it is possible to build the information for other areas of interest. In addition, in order to better integrate the different components (meteorology, chemistry, vegetation) in a fully coupled on-line model, it is preferable to connect the biogenic emissions to the vegetation database already used by the meteorological component.
As BOLCHEM was built merging different codes with different settings policies, and because running the whole model requires the control of a large number of (heterogeneous) parameters that sometimes appear with different names in different parts of the ensemble of setting, the model was equipped with a system that unifies the run setting and allows to manage the whole run, from the definition of the domain (space-time) and resolution to the download of meteo ic/bc, preparation of emissions, and the run itself.
Every parameter is accessible to the user in a unique run.ini with an extensive description that provides a sort user-guide. The BOLCHEM-run system is actually a collection of tools managed by a python script and relies on cfget, a free software specifically developed for the purpose (https://packages.qa.debian.org/c/cfget.html).

Model Setup
The assessment of model performance was performed for a one-year run over Europe (December 2009-November 2010). BOLCHEM was operated in a one-way nesting configuration using two grids on regular lon-lat rotated coordinates where the parent domain extended over Europe with a horizontal resolution of 0.2 • × 0.2 • , while the nested one was centered over Italy at 0.1 • × 0.1 • horizontal resolution (Figure 1). The number of vertical levels is 50 for meteorology and 25 for atmospheric composition (half of the meteorological domain that over flat terrain corresponds to about 5600 m a.s.l) in both simulations. A spin-up period of 30 days (November 2009) was used and not considered for model verification. For the parent run, boundary conditions for the meteorological processor were taken from ECMWF analysis every six hours. Meteorological runs were initialized every day at 00 UTC, while the atmospheric composition was carried through the whole simulation, i.e., the concentration after one-day simulation was used as the initial condition for the next day run. The chemical module (both aerosol and gas phase) are driven by a climatology provided by the Institut national de l'environnement industriel et des risques (INERIS) in the frame of the CityZen project [15].  [8] and Im et al. [9], that evaluated eight operational online-coupled models over Europe and North America for the year 2010 using common emissions and boundary conditions. BOLCHEM results for the Italian domain are compared with the forecast performed by the CAMS regional ensemble model at a horizontal resolution of 0.1 • × 0.1 • . CAMS ensemble model is based on nine (three in 2010) state-of-the-art numerical air quality models developed in Europe. They are combined via an ensemble approach, consisting in calculating the median value of the individual outputs [80].

European Domain
In order to give a general view of the model results, with the purpose to point out the model capability in reproducing the spatial variability over the European and Italian domain, a seasonal average of observed values superimposed to modeled ones over the domain of interest are presented in Figures 2 and 3.    For PM 2.5 and PM 10 the model well reproduces the monthly variability throughout the year with a weak overestimation during spring and underestimation in summer. The larger underestimation in PM 10 than in PM 2.5 , suggests that this could be attributed to Saharan dust storms. Moreover, the model seems to fail to reproduce the PM peaks observed between January and February, probably due to single events of transport at large scale (e.g., Saharan dust) or to accumulation near the ground due to high stability meteorological conditions not well reproduced by the model.
More information on the model capability to describe O 3 behavior, in particular the daily variability, can be obtained by the diurnal cycle, which is mostly driven by solar radiation. O 3 in summer presents a typical diurnal cycle with the maximum located in the early afternoon, which is the main feature affecting air quality. Figure 6 shows, for the summer season, measured concentration vs. modeled concentrations. The model well reproduces the diurnal cycle of ozone concentration, although a slight overestimation during the daytime. Model performance in reproducing the observed concentrations, without any averaging operation (temporal or spatial), is highlighted by the scatter plot diagram. Figure 7 display concentration data as hexagonal points, each having as x-coordinate the measured pollutant concentration (µg m −3 ), and as y-coordinate the simulated pollutant concentration (µg m −3 ) (at the same time and the same Airbase station). In the color bar are reported the number of samples falling in each hexagon. Hourly values have been used for NO 2 and O 3 concentrations, and daily value for PM 2.5 and PM 10 . Results are shown on a seasonal basis. For O 3 , the best agreement is found in the summer. For the spring season, also affected by a moderately high O 3 concentration, the model performances are satisfactory, although a slight model underestimation is evident. More persistent underestimation is present for the cold seasons, especially in winter. This result doesn't represent a serious limit to the model application, being O 3 , unlike PM, not a typical winter pollutant. The concentration of PM 10 is well reproduced by the model in winter, and encouraging results were achieved for the remaining seasons, as pointed out by the basic statistics reported in Table 1.    The simultaneous variation of the correlation coefficient R, of the standard deviation (sigma) and of the (centered) root-mean-square error is shown on seasonal basis by the Taylor diagrams (Figures 8-11).    The diagrams point out that for the O 3 ( Figure 8) the best model performances are in summer and autumn, with the higher value for the correlation coefficient R, the lower value for the RMSE and the normalized SD close to unit. Anyway, also for the other seasons, the value of SD indicates that the variability of the model data is similar to the variability of observations. In winter, the model data having less variability than the measurements (SD < 1) and the correlation coefficient has the lower value. For NO 2 (Figure 9), the model data have more variability than the measurements in all season, except in winter. The correlation coefficient shows similar values in all seasons. Model PM 2.5 data show more variability in summer season ( Figure 10). For PM 10 , Figure 11 highlight the general underestimation of the model data variability in all the season. The model has the best performance in the winter and spring season.
It can however be observed that for all pollutants the model performances are close to those achieved by the current state-of-the-art model system dedicated to air quality study [8,9,80]. Im et al. [8] found that for O 3 the temporal variations were better captured in all models in summer and autumn. The correlation coefficient R was between 0.8 and 0.9 in summer, and between 0.6 and 0.8 in winter. It was found a general underestimation of the annual surface O 3 by up to 18%, especially in winter and spring, while overestimation of the observed concentration was found in autumn. In BOLCHEM the correlation coefficient varies from 0.59 (winter period) to 0.72 (summer period). As regards the mean bias, we found MB = −14.96 µg m −3 in winter and MB = 2.15 µg m −3 in summer. Im et al. [9] showed that for PM the monthly temporal variations were not captured by any of the models, particularly in winter. The rural PM 10 was reproduced with R = 0.18 − 0.86, while the variations at urban sites were reproduced with slightly lower agreement (R = 0.06 − 0.82). Underestimated was found by all models: from 10% to 66% over the rural stations and from 43% to 75% over the urban one. Similar results were obtained for PM 2.5 concentration. In BOLCHEM the correlation coefficient varies between 0.47 and 0.56 for PM 10

Italian Domain
In this section, the results of the nested simulation, with horizontal domain centered over Italy, at resolution of 0.1 • × 0.1 • , are presented and compared with results of CAMS ensemble median. As for the parent simulation, the model capability in reproducing the spatial variability is first considered. Figures 12 and 13 show the seasonal average of observed values at Airbase stations superimposed to modeled ones. Filled circles at Airbase station positions represent the measured concentration. First, we observe that the measurement stations are mainly located in the Po Valley. This area is one of the most polluted regions in Europe, having a high density of anthropogenic emissions combined with frequent stagnant atmospheric conditions. In Figures 12 and 13 the area near big cities, as Milan and Turin, characterized by higher anthropogenic emissions and then higher pollutants concentration value are clearly visible, especially for NO 2 (Figure 12b).  The monthly variability is shown in Figures 14 and 15. As well as for the simulation over Europe, the monthly mean value is considered for all the pollutants, with the exception of O 3 , for which the monthly mean is calculated using the maximum daily 1-h average concentration.  For O 3 (Figure 14a) we found an overestimation during summer, when the O 3 levels are higher. Being the overestimation more evident for the months of July and August (Figure 14a), this could be caused by an overestimation of biogenic VOC emissions, as these are characterized by large uncertainties [81,82]. However, this aspect requires more investigation. The observed overestimation is highlighted by the O 3 daily cycle ( Figure 16). In particular, we can observe a rapid growth during the first hours of the day. This trend could be related to the ability of the model to reproduce the characteristics of the Po Valley (high density of anthropogenic emissions combined with frequent stagnant atmospheric conditions), which makes the greatest contribution to statistics, and not to the higher horizontal resolution, as we can deduce comparing the results obtained with the simulations performed at the two different resolution (Figure 16). The impact of the higher resolution is evident for the night concentration values. Unlike CAMS ensemble median, that overestimates the nighttime O 3 concentration, the nocturnal simulated O 3 concentrations of BOLCHEM are lower than those measured. The O 3 concentration un-derestimation is more pronounced the higher the resolution. One possible cause could be an underestimation of the height of the nocturnal stable boundary layer, with consequent accumulation of NO and therefore removal of O 3 through titration. We recall that the planetary boundary layer is still a hot issue in the study of air pollution and that the estimation of the boundary layer height has an important impact on the dispersion of pollutants. The NO 2 concentration reaches its maximum values in winter, a season in which its average monthly value is well reproduced by the model BOLCHEM (Figure 14b). On the other hand, NO 2 is overestimated in the rest of the year. On the contrary, CAMS ensemble median shows a systematic underestimation in all months. However, we can note that the difference between measured and modeled concentration is smaller, in absolute value, for the BOLCHEM model, as can be seen from the mean bias value ( Table 2). For PM 2.5 ( Figure 15a) and PM 10 (Figure 15b), the model reproduces the seasonal variability, as seen in the data. This variability is not captured by CAMS ensemble median. We also note that in the winter months, more critical due to the high concentration values of the particulate, there is a good agreement between model and measurements.
For providing additional information about the model skill on a seasonal basis, Taylor diagrams are shown in Figures 17 and 18. The basic statistics are reported in Table 2, both for BOLCHEM model and CAMS ensemble median.
O 3 statistics show that BOLCHEM performance is not as good as CAMS ensemble median. However, it can be noted the better performances in summer (with the higher value for the correlation coefficient R) and autumn (with the lower value for the RMSE and the normalized SD close to unit). In case of NO 2 , the best performances are achieved in the winter season in term of MB, even if the correlation is lower than for CAMS ensemble median. As regard PM 2.5 and PM 10 , both in CAMS and BOLCHEM show similar results in summer. For these pollutants, BOLCHEM outperforms CAMS ensemble median in winter, while in the other season the performance is slightly different. (a) (b) Figure 18. As in Figure 17, but for PM 2.5 (a) and PM 10 (b).

Conclusions
We have presented a description of the on-line air-quality model BOLCHEM that is based on the hydrostatic meteorological BOLAM model, the gas chemistry module SAPRC90 and the aerosol dynamic module AERO3. BOLCHEM has been used for different applications, both forecast and air pollution case studies in the framework of different European projects such as GEMS and CESAPO. The model was evaluated on seasonal basis (for the year 2010) for NO 2 , O 3 , PM 2.5 and PM 10 , using measured concentrations available from values measured at background AirBase database. Simulations were performed at 0.2 • × 0.2 • resolution over the European domain, and at 0.1 • × 0.1 • over the Italian domain. The model reproduces the ozone concentration during summer when it has the highest concentrations and impact on vegetation, in particular for the European domain, with correlation coefficient 0.72. The model performances are lower over the Italian domain, for which an overestimation is present during daytime and an underestimation is present during nighttime. The correlation coefficient of O 3 ranges between 0.62 and 0.70, the latter value referred to the autumn period. However, the value of RMSE are comparable with those of some models of the CAMS ensemble, as can be seen from the Evaluation and Quality Assurance reports of these models, provided on a seasonal base (https://atmosphere.copernicus.eu/regional-services). PM 2.5 and PM 10 are particularly well reproduced by the model in wintertime when the PM 10 concentration usually exceeds the limit values imposed by EU Air Quality Directive on ambient air quality and cleaner air for Europe (AQD, Directive 2008/50/EC). Over the European domain for PM 2.5 the correlation coefficient is 0.51 (summer period) and 0.66 (winter period), for PM 10 the correlation coefficient ranges between 0.47 (summer period) and 0.56 (spring period). Comparison of the basic statistics for higher-resolution simulations with those from CAMS ensemble median shows that BOLCHEM better reproduces the aerosol dynamic than CAMS ensemble median. In winter period, for PM 2.5 the correlation coefficient is 0.64 for BOLCHEM and 0.49 for CAMS, for PM 10 the correlation coefficient is 0.63 for BOLCHEM and and 0.56 for CAMS.
Future developments of BOLCHEM will include: 1. estimation of emission potential and foliar biomass density referring to the vegetation dataset used in BOLCHEM, that is UMD Global Land Cover Classification (GLCC); 2.
implementation of a pre-processor system to allow the integration of the mineral dust flux simulated by the model DREAMABOL [83] model into BOLCHEM boundary condition; 3.
upgrade of the gas-phase chemical mechanism and aerosol module; 4.
a major challenge to complete BOLCHEM aerosol feedback description would be to study the effect of aerosol components that are currently not included in the radiation module, such as ammonium and nitrates.