The EMEP MSC-W chemical transport model – technical description

. The Meteorological Synthesizing Centre-West (MSC-W)


D. Simpson et al.: EMEP MSC-W model: description
Convention (now comprising 51 Parties, including USA and Canada) and the European Commission.
The Meteorological Synthesizing Centre-West (MSC-W), based in Oslo, is one of two modelling centres established by EMEP, with main responsibility for photo-oxidant and aerosol modelling. The other modelling centre, MSC-East, is based in Moscow and focuses on heavy metals and persistant organic pollutants. A third Centre, EMEP Chemical Coordinating Centre (CCC) takes care of the EMEP measurement network, and provides the main source of data against which the chemical transport models (CTMs) of EMEP are evaluated (Tørseth et al., 2012). The CTM used at EMEP MSC-W is a 3-D Eulerian model, typically used to tackle problems within the fields of acid deposition, tropospheric ozone, and particles. Results from this model are provided to the International Institute for Applied Systems Analysis (IIASA), providing the atmospheric chemistry results that underpin the GAINS integrated assessment model (http://www.iiasa.ac.at/rains/gains.html).
The last full documentation of the model was Simpson et al. (2003a). Since that time there have been numerous changes, sometimes involving a complete revision of the methodology used (e.g. dry deposition for particles, emissions of hydrocarbons from vegetation, NO emissions from soils, co-deposition of SO 2 and NH 3 , calculation of mixing heights, or introduction of pH response during sulphate formation), and sometimes involving smaller changes in the equations or parameters values. Further, the scope of application of the model has increased enormously. Traditionally, the EMEP model covered all of Europe with a resolution of about 50 km × 50 km, and extending vertically from ground level to the tropopause (100 hPa). The model is now applied on scales ranging from local (ca. 5 km grid size) to global (with 1 degree resolution). The model can now be driven by several different numerical weather prediction or climate models, and has a nesting capability, allowing for example the zooming from 50 km to 5 km scale in the EMEP4UK work of Vieno et al. (2009Vieno et al. ( , 2010. Some of the background for these changes (and some experimental ones) can be found in several recent papers. These include testing of organic aerosol schemes , sea-salt modelling (Tsyro et al., 2011), water-content of aerosols (Tsyro, 2005), ozone deposition Tuovinen and Simpson, 2008), aerosol deposition schemes (Flechard et al., 2011), boundary layer physics (Jeričevič et al., 2010) or soil water modelling (Büker et al., 2011). The chemical schemes mentioned in section Sect. 7 have been compared and a write-up is in progress (Hayman et al., 2012).
The model performance compared to EMEP and other measurements is presented annually in EMEP validation reports, available from www.emep.int, e.g. Gauss et al. (2011). More in-depth discussion can be found in a number of papers. Comparisons for sulphur and nitrogen compounds can be found in Simpson et al. (2006a,b) and Fagerli and Aas (2008). Comparison against trends of inorganic species and EC can be found in Fagerli et al. (2007) and for ozone in Jonson et al. (2006a). Aas et al. (2012) present comparison against AMS and other data-sets during the first so-called EMEP intensive measurement period. Comparisons for large scale CO (and to a lesser extent C 2 H 6 ) have been presented in Angelbratt et al. (2011). The regional forecasts of the EMEP MSC-W model are also constantly under evaluation within the MACC project (Valdebenito and Benedictow, 2011). A discussion of the fine-scale applications and performance of the model can be found in Vieno et al. (2009) and Vieno et al. (2010). Further, the EMEP model has been taking part in a large number of inter-comparisons in recent years (e.g. Cuvelier et al., 2007;Fiore et al., 2009;Huijnen et al., 2010;Jonson et al., 2010a;Colette et al., 2011Colette et al., , 2012Langner et al., 2012).
Given that the EMEP model is being used in a wide range of scientific and policy contexts, there is an urgent need to provide a full description of the model as it is now, and indeed as used in many of the above papers. A short summary of the changes from the 2003 to 2012 model versions can be found in the Supplement, Sect. S1, but the intention of this paper is to present a detailed documentation of the EMEP MSC-W modelling system as it is now. The formulations used by the model are given, along with some details of input data-sets. The aim of this paper is to provide a concise description, rather than discussion, of the model -the latter is left for more extended reports and publications on specific subjects. However, the background to a few of the more recent changes to the model is presented briefly.
Some of the more technical descriptions and tables are provided as a Supplement. For convenience, Table 1 provides an overview of some of the main symbols and abbreviations used in this article. Eliassen et al. (1982) and Eliassen and Saltbones (1983) presented the first long-range transport model within the EMEP framework. The model was Lagrangian, developed for modelling sulphur compounds, and covered the whole of Europe using a 150 km × 150 km grid. This model was further developed for nitrogen compounds (Hov et al., 1988;Iversen, 1990), and ozone (Simpson, 1993(Simpson, , 1995. Eulerian models were subsequently developed for acidification (Berge and Jakobsen, 1998), and photo-oxidants (Jonson et al., 1997(Jonson et al., , 1998(Jonson et al., , 2001. In Simpson et al. (2003a) the first "unified" EMEP model was presented, in which one Eulerian model code was developed for both acidification and photo-oxidant activities.

Short history
In 2008 version rv3.0 of the EMEP model was released as public domain code, along with all required input data for model runs for one year. The second release of the EMEP MSC-W model, denoted EMEP MSC-W version rv3.7 became available in mid 2011, and a new release, rv4.0, is where p * = p S − p T and p, p S and p T are the pressure at level σ , at the surface, and at the top of the model domain (currently 100 hPa), respectively. The model currently uses 20 vertical levels, as illustrated in Fig. 1. The lowest two layers in this system are shown in Fig. 2, with the σ levels from Fig. 1 as solid lines, and the "mid"-layers for which meteorology is generally provided as dashed lines. Diffusion coefficients and vertical velocity, given byσ (= d σ/d t), are valid for the layer boundaries.

The continuity equation
If we let χ represent the mass mixing ratio (kg pollutant per kg air) of any pollutant, the continuity equation may be written: − ∂ ∂σ (σ χp * ) + ∂ ∂σ K σ ∂ ∂σ (χp * ) + p * ρ S The first three terms on the right hand side represent a flux divergence formulation of the advective transport. u, v are the horizontal wind components, and m x , m y are the map factors in the x and y directions (m x = m y in a conformal projection like polar-stereographic). The vertical velocity,σ equals dσ /dt. The 4th term on the right hand side of Eq. (3) represents the vertical eddy diffusion K σ coefficient in σ -coordinates. Horizontal eddy diffusion is not included in the model. In the 5th term, S includes the chemical and other (convection, deposition etc.) source and sink terms.
The numerical solution of the advection terms is based upon the scheme of Bott (1989a,b). The fourth order scheme is utilized in the horizontal directions. In the vertical direction a second order version applicable to variable grid distances is employed.
In our scheme the "air" (χ air =1 kg kg −1 ) is also advected. After each advection step the new mixing ratios are  Fig. 1) and the "mid"-layers for which meteorology is generally provided.
found by dividing the result by the new "air concentrations": (χ x ) t+ t = (χ x p * ) t+ t χ t+ t air , where (χ x p * ) t+ t is the result obtained with the Bott-scheme for component x after a timestep t. This method ensures that, starting with a constant mixing ratio, the result will also be a constant mixing ratio, regardless of the complexity of the wind-fields. The EMEP model's advection scheme is not monotonic, because a monotonicity filter may increase the numerical diffusion. However the scheme will exclude possible negative values of the mixing ratios. The time steps are adapted to the choice of the grid resolution and meteorological data. This work is described in more detail in Wind et al. (2002), and a brief outline is presented in the Supplement, Sect. S2.2.

Convection
An optional (see below) convective mass flux scheme has been implemented in the EMEP model, based on Tiedtke (1989). The implementation is virtually identical to the method used in the Oslo CTM2 model (Berglen et al., 2004), and was originally developed by M. Prather and B. Hannegan, University of California at Irvine (UCI). From the meteorological input data, convective updraft mass flux is provided at every level in each model column and the convective transport of pollutants mass is calculated by the so called elevator principle. The entrainment of air to the updraft cloud core from the surrounding air is calculated as the difference in convective mass flux through the upper and lower boundary of a given grid box, and may be visualised as an elevator stopping at each model layer for air, humidity and pollutant mass to get on or off as illustrated in Fig. 3 (negative entrainment is referred to as detrainment). Vertical transport through convection is much faster than through large scale advection.
As illustrated in Fig. 3 the updraft core will typically gain momentum in the lowest model layers, resulting in a net entrainment, visualised by the upward pointing errors to the left in the lower part of the figure, and lose momentum higher up, resulting in net detrainment. The downdraft core is treated in an analogous way. Within one grid column the downdraft flux is typically about a factor of 10 smaller than the upward flux. The net difference between updraft and downdraft fluxes is treated as a slow subsiding motion. The numerical implementation of the convective routines is described in the Supplement, Sect. S2.1.
Convection is an important process in atmospheric dynamics, but very difficult to parameterise in CTMs (Stevenson et al., 2006). Willett et al. (2008), Zhao et al. (2009) and Monks et al. (2009) (and references cited therein) give examples where significant differences in precipitation and mass transport were found between different parameterisations of convection in NWP models. Used with global-scale IFS meteorological data, the convection module seems to give more realistic results compared to measurements. However, we find that if used with European-scale simulations, we obtain somewhat worse model results compared to observations. This is of course an unsatisfactory situation, but given that all cumulus schemes in NWP models have major uncertainties, we adopt a pragmatic approach and by default switch convection off for the European scale, and on for global scale. The option to switch this module on and off in any case affords some valuable information on the impor-tance of convection, and the uncertainties associated with its implementation.

Nesting
The EMEP MSC-W model now supports 1-way nesting, in which the results of larger-scale runs of the EMEP (or indeed of any other comparable CTM) model can be used as boundary conditions for smaller scale runs. This procedure is most heavily used in the EMEP4UK project (e.g. Vieno et al., 2010), where model runs with 5 km grids over the United Kingdom are nested within larger domain runs of 10 km, which in turn are nested within European scale runs using 50 km grids. (Other configurations are also used). Of course, appropriate meteorological and other data are required for all nesting levels, and for EMEP4UK the WRF (Skamarock and Klemp, 2008) model is used to obtain the necessary data.

Meteorology
During the last few years the EMEP model has been adapted to run with meteorological fields from a number of numerical weather prediction models (NWPs), including PARLAM-PS (Lenschow and Tsyro, 2000;Bjørge and Skålin, 1995;Benedictow, 2003), HIRLAM version version 7.1.3 (Undén et al. 2002, http://hirlam.org/) and ECMWF-IFS Cycle36r1 (http: //www.ecmwf.int/research/ifsdocs/). In 2009 the ECMWF-IFS became available to run with the T799 0.22 • × 0.22 • horizontal spectral resolution and 60 vertical levels on a global domain, and from 2011 we have adopted this model as the default meteorological driver.
For higher resolution modelling, both the EMEP4UK and EMEP4HR projects make use of EMEP model's nesting capabilities (Sect. 2.4) together with the WRF and Aladin models as meteorological drivers -see Vieno et al. (2010), Jeričevič et al. (2010), and associated references for more details.
Meteorological data are normally required at 3-hourly intervals for the EMEP model. Given the wide range of meteorological drivers, which do not all provide all desired model inputs, the EMEP model has systems for deriving parameters when missing, or can do without some meteorological fields. Table 2 summarises the meteorological fields currently used in the EMEP model, and indicates optional fields (one of these, soil moisture index, is briefly discussed in Sect. 3.3). Most 3-D fields are provided at the centre of each model layer, as illustrated in Fig. 2. The horizontal wind components (u and v), and the vertical wind componentσ , are given  Arakawa C-grid (Arakawa and Lamb, 1977). The vertical velocity, given byσ , is provided at the layer boundaries. All other variables are given in the centre of the grid cells. If the vertical wind velocity is not directly available, it is derived from the horizontal wind components and the continuity equations.
Linear interpolation between the 3-hourly values is used to calculate values of these parameters at each advection step. A number of other parameters are derived from these, for example air density, ρ, and the stability parameters and boundary layer heights described below.
Solar radiation is also calculated at every time-step for the deposition calculations, and for photolysis rates, based upon instantaneous values of the solar zenith angle and the model's cloud cover, see Sect. 4.

Boundary layer height (Z PBL )
In general, we characterise the thermal stability of the atmosphere by the bulk Richardson number, which is defined for the layer between any two model levels at heights z j and z k as where g is the acceleration due to gravity, is the difference in horizontal velocity vectors. Following Jeričevič et al. (2010), the mixing height calculation uses a slightly modified bulk Richardson number, Ri B,j , in which z k is always the lowest level (ca. 45 m, cf. k 20 in Fig. 2), but the wind-velocity gradient is referred to ground-level (where V H (0) = 0), thus V H,j,0 = V H,j . The mixing height is defined as the lowest height z PBL = z j at which Ri B,j > 0.25. This formulation is significantly simpler than that used in previous EMEP model versions, and has been shown to provide results which are at least as good (Jeričevič et al., 2010). The method is also very similar to the bulk Richardson number approach used in Seibert et al. (2000), which compared favourably with other methods.
Finally, the PBL height is smoothed with a second order Shapiro filter in space (Shapiro, 1970). The PBL height is not allowed to be less than 100 m or exceed 3000 m.

Eddy diffusion coefficients (K z )
The initial calculation of the vertical exchange coefficients (K z , units m 2 s −1 ) is done for the whole 3-D domain, using: where the critical Richardson number Ri c is given by: Ri crit,k = A ( z k / z 0 ) B , A=0.115, B=0.175 and z 0 =0.01 m (Pielke, 2002), is the turbulent mixing length, and V H represents the difference in wind-speed between two grid-cell centres separated by distance z, and K min = 0.001 m 2 s −1 . The numerical values follow from the suggestions of Blackadar (1979) and Pielke (2002).
The turbulent mixing length, , is parameterized according to: = k z , z ≤ z m = k z m , z > z m where k is the von Karman's constant (0.41), z is the height above the ground and z m = 200 m.
Below the mixing height z PBL , these K z values are recalculated. For neutral and stable conditions the simple formulation of Jeričevič et al. (2010) is used, whereby: for z < z PBL . The values 0.39 and 0.21 are empirical constants derived from large eddy simulation experiments. u * is the friction velocity provided by the NWP model (= √ τ/ρ, m s −1 ).
For unstable situations, new K z values are calculated for layers below the mixing height using the O'Brien (1970) profile: in which h s is the height of the surface layer (or the socalled constant flux layer), which we set to be 4 % of the mixing height z PBL (Pielke, 2002). From the similarity theory of Monin-Obukhov (e.g. Stull, 1988;Garratt, 1992) we have where h is the atmospheric stability function for temperature, assumed valid for all scalars. The latter is derived using standard similarity theory profiles (Garratt, 1992). The Obukhov length is given by: where c p is the specific heat capacity of dry air (1005 J kg −1 K −1 ), and ρ is air density. The sign here is consistent with H directed away from the surface (positive H gives unstable conditions). Finally, in sigma coordinates, the diffusion coefficient has the following form:

Soil water
Soil water (SW) is very difficult to model accurately in largescale models, since it depends very much on assumptions concerning parameters such as soil texture, and vegetation characteristics such as rooting depth that are not generally amenable to measurements (e.g. Baker, 2003;Büker et al., 2011;Miller et al., 2007). Different NWP models also make use of very different schemes for soil water, depending on the complexity of the underlying vegetation schemes, and these models provide different outputs -sometimes SW in terms of volumetric amount (e.g. from HIRLAM), sometimes in terms of a soil moisture index (ECMWF, discussed below). Volumetric outputs can be difficult to interpret unless the associated soil and vegetation characteristics are known for that NWP. Soil moisture is important though for dry-deposition and dust emission rates, so we have implemented a procedure which unifies the treatment from different NWP models. The exact methodology depends on the NWP model and its SW outputs, but essentially we define minimum and maximum soil water amounts to be SW min (identified with wilting points for example) and SW max (identified with field capacity), which may be constant over an NWP domain, or vary spatially, and then define the soil moisture index (which we previously denoted as relative extractable water index), as: The index S MI has the advantage over volumetric methods that it is less sensitive to local soil characteristics, and hence is easier to interpolate across different vegetation types and grids. For example, a reasonable estimate of volumetric SW can be made given local values for SW min and SW max , if S MI is known. The ECMWF IFS data we now use by default provides S MI values directly; these are available for the near-surface (ca. 10 cm) and deeper (1 m) soil layers, which we use for dust and dry-deposition modules, respectively.

Radiation
Calculation of direct and diffuse radiation is needed for chemical photolysis rates (Sect. 7.3), and in addition, calculation of photosynthetically active radiation (PAR) is needed for calculating biogenic VOC emissions (Sect. 6.6), and for calculation of stomatal conductance for dry deposition or ozone uptake modelling (Sect. 8).  For radiation calculations at level k in the model, we need an estimate of the integrated cloud fraction in the column above and including k. We use a maximum overlap assumption, in which the fraction f k cloud is set to the maximum value of the cell-volume cloud covers from k and all layers above, i.e. from 1...k, cf. Fig. 2.
Following Pierce and Waldruff (1991) and Iqbal (1983), direct normal irradiance (W m −2 ) is then estimated as: where C N is a clearness number, assumed equal to 1, T k is a transmissivity factor (set as T k = 1 − 0.75f 3.4 cloud ), A, B are empirical co-efficients from Iqbal (1983), θ is the solar zenith angle, p is the local pressure (Pa) and p 0 is standard sea-level pressure, set equal to 101.3 kPa.
The direct and diffusive radiation on a horizontal surface (W m −2 ) are then given simply by: where the co-efficient C is also taken from Iqbal (1983).
Calculation of PAR values are made for each vegetated land-cover class within the grid, as PAR depends on the canopy's leaf area index (LAI). Following Norman (1979Norman ( , 1982 we divide the canopy into sunlit and shaded leaves, and calculate the leaf-area and PAR for each class with: LAI shade = LAI − LAI sun (15) where α is the average inclination of leaves in the canopy (assumed 60 • to represent a spherical leaf distribution).

Land-cover
Land-cover data are required in the model, primarily for dry deposition modelling and for estimation of biogenic Atmos. Chem. Phys., 12, 7825-7865, 2012 www.atmos-chem-phys.net/12/7825/2012/ emissions. As noted in Sect. 2, the standard EMEP grid has a resolution of approx. 50 km × 50 km, but grid sizes in reported applications have ranged from 5 km × 5 km to 1 • × 1 • . Whatever the size, the land-use databases give the fractional coverage of different land-cover types within each surface grid cell. This allows sub-grid modelling using a socalled mosaic approach -allowing for example ecosystem specific deposition estimates.
The 16 basic land-cover classes are summarised in Table 3. Additional land-use classes are easily defined and indeed the specific categories "IAM DF", "IAM MF" and "IAM CR" are assigned for provision of data to ozone-effects studies and integrated assessment studies (e.g. Mills et al., 2011a,b). For European scale modelling the land-cover data are derived from the CORINE system and from the Stockholm Environment Institute at York (SEIY) system (www.york. ac.uk/http://www.sei-international.org/landcover). The basic principle used was to apply CORINE data wherever available, thereafter SEIY data. In addition, the more detailed SEIY data (especially on agriculture) was used to guide the split of the broader CORINE categories into the EMEP landclasses needed by the model. The final merge of these data was done at the the LRTAP Coordination Centre for Effects (CCE at RIVM, Posch et al. 2005). For global scale runs, land-cover from GLC-2000 (http://bioval.jrc.ec.europa. eu/products/glc2000/glc2000.php) are used.
For the vegetative land-cover categories for which stomatal modelling is undertaken (see Sect. 8.5), a number of phenological characteristics are needed. By default, these are specified in input tables for each EMEP land-cover c . In particular, the start and end of the growing season (SGS, EGS) must be specified. The development of leaf area index (LAI) within this growing season is modelled with a simple function as illustrated in Fig. 4. The parameter values used for these LAI estimates are given in Table 3.

Emissions
The standard emissions input required by EMEP model consists of gridded annual national emissions of sulphur dioxide (SO 2 ), nitrogen oxides (NO x = NO + NO 2 ), ammonia (NH 3 ), non-methane volatile organic compounds (NMVOC), carbon monoxide (CO), and particulates (PM 2.5 , and PM c , the latter being the coarse aerosol fraction, PM 10 -PM 2.5 ). The particulate matter categories can be further divided into elemental carbon, organic matter and other compounds as required. Emissions can be from anthropogenic sources (burning of fossil and biomass based fuels, solvent release, etc.), or from natural sources such as foliar VOC emissions or volcanoes. Several sources are hard to categorise as anthropogenic versus natural , e.g. with emissions of NO from microbes in soils being promoted by N-deposition and fertilizer usage. 78 Fig. 4. Schematic of LAI development and associated parameters. SGS and EGS are the start and end of the growing season, in day-numbers. L S and L E represent the length of the LAI-increase and decline periods, also in day-numbers. Maximum and minimum (within the growing season) LAI values are given by LAI max , LAI min .

National EMEP emissions
As part of the EMEP Protocol under CLRTAP, national estimates of the anthropogenic emissions should be provided to EMEP every year from each country, along with spatial distribution to the EMEP grid. These emissions are provided for 10 anthropogenic source-sectors denoted by socalled SNAP codes. An eleventh source-sector exists in the officially-submitted database ("Other sources and sinks"), but this consists almost entirely of emissions from natural and biogenic sources. Officially submitted emissions from such sources are not used in the modelling work, except for those from volcanoes (sections 6.6-6.11 discuss the methods used for dealing with such emissions in the modelling framework). Further details of the anthropogenic emissions can be found in Mareckova et al. (2009); the emission database is available from http://www.emep.int, and further details can be obtained at that site. Figure S1 in Supplement illustrates the spatial distribution of two sets of data for these anthropogenic emissions (NO x and SO 2 ), and two sets of data for biogenic VOC emissions.

Vertical distribution
These land-based gridded emissions are distributed vertically according to a default distribution based upon the SNAP codes, as shown in the Supplement, Table S3. These distributions were originally based upon plume-rise calculations performed for different types of emission source which are thought typical for different emission categories, under a range of stability conditions (Vidic, 2002), but have since been simplified and adjusted to reflect recent findings (Bieser et al., 2011;Pregger and Friedrich, 2009). The biggest change has been in sector 2 (non-industrial combustion), where now 100 % of the emissions are placed in the lowest model layer, reflecting the large dominance of domestic combustion for this emission category.

Temporal distribution
For most SNAP sectors, emissions are distributed temporally according to: where f i t is the temporal factor for SNAP sector i, country c, and f i,c m,y , f i,c d , f i h , are factors accounting for month (and year for SNAP-1, see below), day-of-week (or for SNAP-2 day of year, see below) and hour of the day. These factors are derived largely from data provided by the University of Stuttgart (IER) as part of the GENEMIS project (Friedrich and Reis, 2004), and are available as data files from the EMEP model website, www.emep.int. They are specific to each pollutant (except f i h ), emission sector, and country, and thus reflect the very different climates and hence energy-use patterns in different parts of Europe. Fig. 5 illustrates the monthly variations in emissions of SO x , NO x and NMVOC for selected countries in different parts of Europe. The annual cycles of SO x and NO x emissions for this year (2008) are rather constant in the Western European countries (Sweden, UK, Spain), but still show winter peaks in the two Eastern European countries (Poland, Ukraine). The ratio of SO x to NO x emissions varies markedly from country to country. Note that these plots illustrate total emissions, and the flat cycles for SO x and NOx may partly be ascribed to the importance of traffic emissions (which have very low seasonal cycles), and the lower winter/summer ratio assumed for SNAP-1 in recent years (below). Emissions for particular sectors can show much stronger variation; an example for the domestic emissions of organic carbon emissions can be found in Bergström et al. (2012). The spatial distribution of BVOC emissions is presented in the Supplement, Fig. S1.
The three improvements which have been made to this methodology in 2011-2012 versions of the model are discussed below:

SNAP-1: decreasing winter/summer ratios
The temporal patterns from GENEMIS were derived for the year 1994, and prior to rv3.9 these values were used for all years, i.e. f i,c m was the same set of 12 values for all years. However, as illustrated in Grennfelt and Hov (2005), the winter/summer ratios of electricity consumption have been decreasing in recent years, from about 1.33 in 1990 for the UK to 1.22 in 2000, and from 1.1 to 1.02 over the same period in Italy. Despite very different climates, these changes both represent a 10 % decrease in the winter/summer ratio over these 10 years. Discussion with IIASA (M. Amann personal communication, 2011) suggest that this decreasing trend has continued. For SNAP-1, power stations and suchlike, we therefore modify these variations, "flattening" the monthly factors towards the annual mean by a factor ranging from 0-10 % between 1990 and 2010: where f m,1994 is the monthly factor obtained from GENE-MIS for 1994, (y) is set to zero before 1990, y − 1990 between 1990 and 2010, and to 1 after 2010. The cosine term provides an annual cycle, and m − 8 ensures that maximum changes occur in February and August. (Note that the mean of all f m,y factors is 1.0, we are just changing the amplitude of the annual cycle.)

SNAP-2: use of degree-day factors
SNAP-2 consists mainly of domestic combustion, and as of rv3.9, the day to day variation is based upon a modification of the heating degree day concept. For day of year j , with mean daily temperature in • C of T 24h j we set the heating degree day to be H dd,j = max(18 − T 24h j , 1). (The minimum value of 1 is used rather than zero just to avoid numerical problems). These degree days are pre-calculated in the model for each grid cell, and averaged to find the annual mean degree-days for each grid-cell, H dd .
These degree-day factors are so far country-independent, being a function only of gridded daily temperatures. However, the GENEMIS monthly factors for SNAP-2 are used to establish a minimum 'base' factor for each country, f c B , which in some countries would include summertime use of gas-appliances for cooking, etc. Time-variation of emissions above this level are driven by calculations of heating degreedays. For day-number j , SNAP-2 we have: Thus, in summertime where temperatures are close to or exceed 18 • C, this emission factor is very small, but in winter the factor is usually significant, and can change quite substantially from day to day.

Hourly emissions
Earlier versions (up to rv4β) of the EMEP model used simple day-night factors (see Table S4) to allocate emissions within the day. In version rv4 we make use of new hourly data, provided by B. Bessagnet, INERIS, as part of ongoing work for the EMEP task Force on Measurements and Modelling. These data consist of a matrix of 11 SNAP sectors ×7 days per week × 24 h. These values are somewhat simplified versions of the hourly data presented by Menut et al. (2012).

VOC speciation
Speciation of VOC emissions is also specified separately for each source-sector. The EMEP model uses a "lumped-molecule" approach to VOC emissions and modelling, in which for example the model species n-butane represents all C3+ alkanes, and o-xylene represents all aromatic species (Andersson-Sköld and Simpson, 1997). As discussed in more detail in Hayman et al. (2012), the VOC data used in the current EMEP model are derived from the detailed United Kingdom speciation given in Passant (2002). Although the exact VOC speciation used can be varied to suit particular emission scenarios (e.g. Reis et al., 2000), the default split is typically used, as given in the Supplement, Table S5.

PM speciation
Where elemental and/or organic carbon (EC, OC) are required, emissions of PM 2.5 and PM 10 need to be speciated into these components. In fact, we are often interested in emissions of organic matter, OM, which includes for example oxygen, hydrogen and other atoms bound to the OC. In order to generate these speciations, we make use of country specific information on EC, OC and PM emissions provided by IIASA. For the fine PM fraction, OM emissions by mass are assumed to be 1.3 times the OC emission, although with a cap to make sure that EC + OM ≤ 0.99 PM. For the even more uncertain coarse fraction, we use a simple default speciation as given in the Supplement, Table S6. For some studies, explicit emissions of EC (or related black carbon, BC) have been available, e.g. for the modelling studies within the CARBOSOL project (Fagerli et al., 2007;Simpson et al., 2007b;Tsyro et al., 2007) emissions from Kupiainen and Klimont (2007) were used, and for the EU-CAARI project (Kulmala et al., 2011;Bergström et al., 2012) emissions were from van der Gon et al. (2009).

Aircraft
Emissions of NO x from aircraft are provided by data from the EU-Framework Programme 6 Integrated Project QUAN-TIFY. The data have been downloaded from the project website www.pa.op.dlr.de/quantify. The emissions are calculated on an annual basis and disaggregated according to a seasonal variation to create monthly files on a spatial resolution of 1 • × 1 • × 610 m. The emissions are interpolated to the relevant model grid during model runtime. In the EMEP model, only NO x emissions from aircraft are used so far.

Shipping
The emissions from international shipping were created originally by ENTEC (now part of AMEC Environment Infrastructure, UK, www.amec-ukenvironment.com) and IIASA, and recently in the context of the revision of national emission ceilings directive as described in Cofala et al. (2007) and Jonson et al. (2009). The latest data take account of reduced sulphur emissions in recent years. Data are now available for NO x , SO x and PM (for 2000, 2005, 2010, 2015, 2020, 2025 and 2030), with interpolation between these years when required.
Emissions from national shipping are not included in this ship inventory as national emissions should be included in the reported emissions (SNAP sector 8) to UN-ECE by the individual parties to LRTAP Convention. Unfortunately not all countries report emissions from national shipping, and for those who do it can not be distinguished from other mobile sources.

Foliar NMVOC emissions
Biogenic emissions of isoprene and (if required) monoterpenes are calculated in the model for every grid-cell, and at every model time-step, using near-surface air temperature (T 2 ) and photosynthetically active radiation (PAR, see Sect. 4). Following the ideas proposed in Guenther et al. (1993Guenther et al. ( , 1995, the first step in the emission processing is to define "standard" emission factors, which give the emissions of particular land-covers at standard environmental conditions (30 • C and PAR of 1000 µmol m −2 s −1 ).
Emission factors for forests have been created from the the map of forest species generated by Köble and Seufert (2001). This work (also used by Karl et al. 2009 andKesik et al. 2005) provided maps for 115 tree species in 30 European countries, based upon a compilation of data from the ICP-forest network UN-ECE (1998). These data were further processed to the EMEP grid (S. Cinderby, SEIY, personal communication, 2004).
The EMEP model cannot deal with all these different forest species, but rather has maps of aggregated land-cover types, such as temperate/boreal coniferous forest (CF), as in Table 3. Emission rates for the EMEP aggregated landcover classes ( c ) are developed from maps of the Köble and Seufert (2001) land-cover types (λ c ) with: where E * c ,i is the area-specific reference emission rate (µg m −2 h −1 ) for an EMEP land-cover class, at standard environmental conditions, ε * λ c,i is the mass-specific emission rate (µg g −1 (dry-weight) h −1 ) for BVOC compound i and a particular real land-cover class (λ c ) at these standard conditions, A λ c is the area, and D λ c is the foliar biomass density of that species. The delta (δ) function is set to 1.0 for those species (λ c ) belonging to the EMEP land-cover group ( c ), zero otherwise. The standard emission factors are as given in the Supplement , Table S7.
For example, the standard emissions factor for the CF landcover (temperate/coniferous forest) would be calculated as the weighted sum of the species-specific emissions factors for any species included in this category, thus c would include Norway spruce, Sitka spruce, Scots pine, etc. The 57 γ L and γ T,iso as in Guenther et al. (1993) MTP Guenther et al. (1993) MTL =γ L,iso =γ T,MTP 0.57 Light-dependent monoterpene emissions Notes: all coefficients from Guenther et al. (1993), resulting E * c,i give standard emission factors per m 2 of the appropriate EMEP landcover category.
These E * c,i maps are intended to represent broad species characteristics rather than to capture details of the spatial distribution, and in order to reflect this we have smoothed the emission factor fields using a simple distance weighted filter.
For non-forest vegetation types (e.g. grasslands, seminatural vegetation) or for forest areas not covered by the emission factor maps described above (e.g. for eastern Russia, or non-European forests when modelling at global scale), default emission factors are applied. These factors are given in Table 3.
Emission potentials are then re-calculated to instantaneous emissions every time-step in the model (every 20 min), using the grid-cell relevant temperature and radiation conditions: where E c ,i is the temperature and (where appropriate light) corrected emission per square meter of EMEP landcover c . The environmental correction factor γ c ,i consists of corrections for the canopy LAI, temperature, light and canopyshading: where the LAI factor, γ LAI is simply defined as LAI/LAImax for each land-cover c . In the EMEP model, γ CAN,i accounts for the effects of shading throughout the canopy. In principle a multi-layer canopy model could be used to specify leaf temperature and radiation conditions at different vertical levels. However, here we use a simple non-canopy approach, assuming that ambient temperature is similar to leaf temperature and that the use of "branch-level" emission potentials, which are typically a factor 1.75 smaller than leaf-level values (Guenther et al., 1994), accounts for the shading effect. Tests in European conditions have suggested differences in total emis-sions between the two methodologies of around 20 % . Given the many uncertainties introduced by the forest-canopy model itself (e.g. in temperature and light profiles within the canopy), and the lack of evaluation of such models under European conditions, we use the same procedure as Simpson et al. (1999) and simply specify that γ CAN,i = 1/1.75 = 0.57 for light-sensitive emissions and γ CAN,i = 1 for the pool terpenes.
The light correction factor γ L and temperature correction factor γ T are different for the model's three emission categories: isoprene, pool-dependent monoterpenes (MTP) and light-dependent monoterpenes (MTL). Isoprene is always light and temperature controlled. MTP emissions are derived entirely from pool-emissions, and so have γ L = 1 always. MTL emissions are synthesised, and are both light and temperate controlled. Table 4 summarizes the environmental correction factors used. Figure 5 illustrates the monthly variations in emissions of isoprene and monoterpenes for selected countries in different parts of Europe, also in comparison to the anthropogenic emissions and soil-NO (below). These results clearly illustrate not just the strong seasonal cycle, but also the large country to country differences. In the UK for example, BVOC emissions are smaller than anthropogenic even in the summer months, but in the other countries summertime BVOC emissions can be far greater than anthropogenic NMVOC. Emissions of monoterpenes dominate over those of isoprene, also in most of the countries that are not shown. Annual emissions of these BVOC are given in Table S2 of the Supplement, where again the importance of these sources is obvious. The spatial distribution of isoprene and monoterpene emissions are shown in Fig. S1  For global scale modelling the EMEP model can make use of monthly averaged soil NO emissions from a process-based terrestrial-biosphere model (Zaehle et al., 2011), kindly provided as netcdf files with 1 • × 1 • resolution (S. Zaehle, personal communication, 2010).
For European-scale applications, we make use of more detailed land-cover and meteorological data. Emissions of NO from soils of seminatural ecosystems are specified as a function of the N-deposition and temperature: where E * NO, c is the maximum emission rate, set to 150 µg(N) m −2 h −1 for coniferous forest, and 50 µg(N) m −2 h −1 for deciduous forests and other seminatural ecosystems. N T is the temperature response, identical to that used by Laville et al. (2005) and Linn and Doran (1984), and which also seems broadly consistent with data presented by Schaufler et al. (2010). f N dep is a scaling factor to account for the N-deposition load in each grid. For f N dep we take the ratio of annual deposition divided by 5000 mg(N) m −2 , with maximum value 1.0.
For crops, emissions are given by: where E * NO, c is 80 µg(N) m −2 h −1 for all crops, The function f β,n d applies a β(2, 2) function, which produces a value 1.0 when the daynumber n d (between 1 to 366) is equal to the start of the growing season (SGS), falling to zero 30 days on either side of SGS. E 0 NO is the baseline emission level of 1 µg(N) m −2 h −1 .
The approaches used are meant to loosely capture two of the most important dependencies found in field and experimental studies. For example, from a detailed study of 15 forest sites across Europe, Pilegaard et al. (2005) found an almost linear relationship between NO emissions and Ndeposition at coniferous sites, with emissions ranging from non-detectable at a Finnish site to ca. 80 µg(N) m −2 h −1 at two high-deposition sites in the Netherlands and Germany. For deciduous forests the relationship with N-deposition was much weaker, with rates varying from 0.7 (Scotland) to 13 µg(N) m −2 h −1 (Germany). The deposition estimates were based upon throughfall for coniferous forest, and throughfall plus stem-flow for deciduous, and so are both uncertain and not strictly comparable. Schaufler et al. (2010) found a somewhat closer relationships between soils from coniferous and deciduous forests in an experimental study, albeit with only a few sites.
The procedure used for crops is designed to loosely mimic results shown in for example Butterbach-Bahl et al. (2009), Rolland et al. (2008Rolland et al. ( , 2010, or Laville et al. (2005Laville et al. ( , 2009), all showing a broad peak in emissions in springtime (corre-sponding to the application of fertilizer and start of the growing season). Figure 5 illustrates the monthly variations in emissions of soil-NO in comparison to the anthropogenic emissions of NO x . The spring peak is clearly seen, starting earlier in southern compared to northern Europe. Country to country differences are large. For Sweden for example, a heavily forested country with relatively low population density, soil-NO emissions are rather large compared to the (low) anthropogenic emissions. For the densely populated United Kingdom, on the other hand, the soil-NO emissions are almost negligible compared to those from industry and traffic. Annual emissions of these soil-NO are given in Table S2 of the Supplement. It should be noted that, although relatively small in most countries (compared to the combustion-sources of NO x ), these emissions can still impact atmospheric chemistry because of their seasonal cycle and location in NO xsensitive areas Butterbach-Bahl et al., 2009).
This methodology has of course a number of weaknesses, including lack of controls by soil moisture, but the emission rates seem to correspond reasonably well to the (widelyscattered) published values from European forests and agricultural areas cited above. A more detailed methodology would require data on a host of factors which are not normally available at the European scale, including details of soil and vegetation types, and timing of crop growing seasons, fertilization, and irrigation.

Sea salt
The generation of sea salt aerosol over oceans is driven by the surface wind. There are two main mechanisms for sea salt aerosol generation: bubble bursting during whitecap formation (indirect) and through spume drops under the wave breaking (direct). The latter mechanism is believed to be important source for particles larger than 10 µm and at wind speeds exceeding 10-12 m s −1 . In the EMEP MSC-W model, sea salt calculations include primarily particles with ambient diameters up to 10 µm. These sea salt particles originate mainly from the bubble-mediated sea spray. As discussed in detail in Tsyro et al. (2011), the EMEP model's parameterisation scheme for calculating sea salt generation is based on two source functions, those of Monahan et al. (1986) and Mårtensson et al. (2003). The equations used are briefly recapitulated in the Supplement (Sect. S4.5), but the reader is referred to Tsyro et al. (2011) for a thorough discussion and comparison with measurements and other models. resolution, on a fine 1 km × 1 km grid. We store these data on a slightly coarser grid (0.2 • × 0.2 • ) globally for access by the EMEP model.
For earlier years, and in previous versions of the model (e.g. as used in Hodnebrog et al. 2012, the model used the 8-daily fire emissions from GFED-2 (Global Forest Emission database, http://www. globalfiredata.org), as documented in van der Werf et al. (2010).
Emissions from either database include SO 2 , CO, NO x , NMHC, PM 2.5 , PM 10 , OC, and BC. Where OM is needed explicitly, we scale from OC using a factor of 1.7 (based on AMS measurements presented by Aiken et al. 2008). Emissions are homogeneously distributed over the eight lowest model layers, loosely following recommendations by Sofiev et al. (2009) to use a PBL height as an approximate height for emission injection.

Dust
The sources of natural mineral dust in the model include windblown dust from deserts, semi-arid areas, agricultural and bare lands within the model domain, as well as dust produced beyond the model grid (e.g. on African deserts) and transported to the calculation domain. A preliminary road-dust module has also been implemented. The inflow of African dust is accounted for through boundary conditions. The monthly average concentrations of fine and coarse dust, calculated with the global chemical transport model of the University of Oslo (CTM-2) for 2000, are currently used as boundary conditions (Grini et al., 2005).
The parameterisation of wind mobilisation of soil dust is based among others on the works of Marticorena and Bergametti (1995), Marticorena et al. (1997), Alfaro and Gomes (2001), Gomes et al. (2003), and Zender et al. (2003). The key parameter driving dust emissions is wind friction velocity. The dust mobilisation by wind and the horizontal motion of soil particles (called saltation) occurs when the wind friction velocity exceeds a threshold value. This threshold value depends on the size of soil aggregates. The model employs a partitioning scheme of wind shear stress between the erodible and non-erodible surface elements to calculate the threshold friction velocity (Marticorena and Bergametti, 1995). Currently, the threshold friction velocity is calculated for a particle size optimal for saltation, which is assumed to be 75 µm (Zender et al., 2003). The general expression for threshold wind friction velocity (u * ,th ) is written as where u * ,sm is the threshold friction velocity for erodible (smooth) part of surface, f eff is the efficient friction velocity ratio (describing wind drag partitioning between erodible surface and non-erodible roughness elements), and f w is the correction factor accounting for soil moisture. Following Marticorena and Bergametti (1995), where z 0,s is the roughness length of the erodible part of the surface (smooth), i.e. roughness of soil aggregates, z 0 is the roughness length of the non-erodible roughness elements (e.g. pebbles, rocks, vegetation). The roughness length of smooth erodible surfaces depends on soil morphology and is calculated following Marticorena et al. (1997) as z 0,s = d s /30, where d s is the diameter of erodible particles, for which the median diameter of the most coarse population of the soil is used. The suppression of soil erosion by soil moisture is accounted for as suggested by Fécan et al. (1998). The correction factor accounting for increase of threshold friction velocity due to soil moisture is calculated as where w is the gravimetric soil moisture (kg kg −1 ) and w is the minimum soil moisture from which the threshold velocity increases. The latter depends on soil texture as: where F clay is the fractional clay content of the soil. In the present model version, w is assumed to be equal to the Permanent Wilting Point (PWP) obtained from ECMWF-IFS data. The EMEP model's soil moisture index (Sect. 3.3) is converted first to volumetric, and then gravimetric soil water, following Zender et al. (2003), using information on sand content in the soil. The land-use types, from which windblown dust emissions are calculated in the model, include deserts/bare lands and agricultural arable lands outside growing periods. Some additional constraints are imposed on the onset of windblown dust generation, so that no emissions take place: (1) during precipitation events (with precipitation rate greater than 0.2 mm per day) and two days afterwards; (2) under high surface relative humidity conditions (RH > 85 %); and (3) from frozen surface or surface covered by snow.
The condition for the onset of dust mobilisation by wind is u * ≥ u * ,th . The model allows a possibility of accounting for the gustiness of wind at free convection conditions. As proposed in Beljaars (1994) where V 10 is the velocity of horizontal wind at 10 m height, w * is the free convection velocity scale, z 0 is the land-use defined roughness length and z 10 = 10 m. The term 1.2w 2 * represents the near surface wind induced by large eddies. The horizontal flux of soil particles (i.e. saltation) is calculated as in Marticorena and Bergametti (1995) where Q s is the horizontal mass flux of soil particles (kg m −1 s −1 ), ρ air is the air density, g is the gravitational acceleration and C is the empirical coefficient (C = 2.61 based on Zender et al. (2003) and references therein). The vertical flux of dust particles, released by sandblasting mechanism from the saltating and/or surface soil aggregates, is simulated as where F is the vertical mass flux of dust (kg m −2 s −1 ), A s is the area fraction of erodible soil in the grid cell, K is the coefficient accounting for soil erodibility (or availability of loose soil aggregates), and α is the sandblasting efficiency (m −1 ). Based on the experimental results in Gomes et al. (2003), the following values (providing the best fit with measurements) are currently used in the model: α = 2.0 × 10 −5 , 1.5 × 10 −5 and 1.0 × 10 −5 m −1 and K = 0.5, 0.05 and 0.02 for North African deserts, Mediterranean arid areas and arable lands respectively. Finally, we can mention that a road dust emission module has been added to the EMEP model. The code and methodology are taken directly from that described in Denier van der Gon et al. (2010). This is very preliminary work, however, and details of the results will be presented elsewhere.

Other sources
Biogenic emissions of dimethly sulphide (DMS) can sometimes make a significant input to European sulphate levels. As discussed in detail by Tarrasón et al. (1995), the EMEP model uses a very simplified treatment, in which DMS is not modelled explicitly, but rather we assume that most DMS enters long-range transport already as sulphur dioxide. Monthly emission fields of DMS-derived SO 2 are taken from the work of Tarrasón et al. (1995).
Emissions of volcanoes are introduced into the model as point sources, at a height determined by the height of each volcano. For the standard European-scale runs, volcano emissions are based upon officially reported data. These have been provided by Italy for many years, and recently by Iceland. For global and regional scale calculations, a new module for volcanic eruptions with default values based upon Mastin et al. (2009a,b) has recently been implemented and is currently in testing.
Emissions of NO x from lightning are included as monthly averages of global 3-D fields on a T21 (5.65 • × 5.65 • ) resolution (Köhler et al., 1995).

Chemistry
The chemical scheme used for gas-phase chemistry traces its origins to the EMEP chemical mechanisms that began with Eliassen et al. (1982). This scheme has been updated and tested against other schemes in a number of studies Simpson, 1995;Kuhn et al., 1998;Andersson-Sköld and Simpson, 1999). The scheme documented in Simpson et al. (2003a) and Andersson-Sköld and  is now denoted EmChem03. The latest scheme was largely developed in 2008-2009 and is denoted EmChem09. Compared to EmChem03, EmChem09 has updated rate-coefficients, and some additional species, including HONO. A detailed comparison of these chemical schemes, including their response to emission changes is presented in Hayman et al. (2012).
The EMEP model now uses a chemical pre-processor (GenChem) to convert lists of input chemical species and reactions to differential equations in Fortran code. At the time of writing, eight different chemical schemes have been tested within the EMEP model, as discussed in detail in Hayman et al. (2012) and summarised in Table 5. A large number of schemes for organic aerosol have also been tested (Simpson et al., 2007b;Bergström et al., 2012), but these are too complex and numerous to document here. Here we document only the default chemical scheme, EmChem09.
The numerical solution of the chemical equations is discussed in Sect. 7.8 and Supplement, Sect. S2.3.

Species used, EmChem09
Tables S8-S10 list the chemical compounds used in the Em-Chem09 scheme, along with associated characteristics such as the assignments used for dry and wet deposition. Most species are sufficiently long lived that they are included in both the advection and chemical equations. The species labelled "short-lived" have sufficiently short lifetimes that their concentrations are essentially controlled by local chemistry, so they are not included among the advected species. (Some short-lived species are advected anyway for numerical reasons.) Note that this list excludes a number of intermediate species which are assumed to react immediately upon Atmos. Chem. Phys., 12, 7825-7865 formation. For example, H atoms react immediately with O 2 to form HO 2 , and so are not included explicitly. The EMEP model distinguishes five classes of fine and coarse particles, which for dry-deposition purposes are assigned mass-median diameters (D p ), geometric standard deviations (σ g ), and densities (ρ p ). The characteristics of these aerosol classes are given in Table 6.
It can be noted that the assumed D p for coarse nitrate particles has been reduced in recent years compared to Simpson et al. (2003a) which had D p = 4 µm. This choice reflects an assumption that coarse nitrate formation is driven by surface-area rather than mass (hence favouring the smaller size-ranges), and consistent with Pakkanen et al. (1996) and Torseth et al. (2000). This assumption is very uncertain however, and probably depends on whether dust or sea-salt is the main reacting surface. In future we will consider explicit modelling of nitrate formation on different types of aerosol in order to better characterise the size-distribution. Accounting for the difference betweem MMD and aerodynamic diameter, the choice that D p = 3 µm implies that 27 % of calculated coarse-nitrate can be assigned to the PM 2.5 fraction.
The semi-volatile organic compounds involved in SOA formation are a special case, in that the model transports both the gas and the aerosol fraction as one lumped concentration for numerical stability. The model also tracks the gas fraction as a separate quantity. For these compounds, dry and wet deposition processes are applied as appropriate to the different fractions.

Gas-phase chemical mechanism
Table S11 lists the chemical reaction mechanism used in the photo-oxidant model (for photolysis reactions, see below). Rate-coefficients for 3-body and some other reactions are given in Tables S12-S13. During 2008-2009 the scheme's rate-coefficients have been updated and in some cases replaced by Troe expressions to allow their application to the greater range of temperatures and pressures inherent in the 3-D model domain. The rates and products were updated to be, as far as possible, consistent with IUPAC recommendations (http://www.iupac-kinetic.ch.cam.ac.uk/); most of the reaction coefficients are from Atkinson et al. (2004Atkinson et al. ( , 2006. Table S14 lists the photolysis reactions used in the model for the EmChem09 mechanism. The reactions are taken from Simpson et al. (1993), with minor updates. The calculation of photo-dissociation rates (J-values) is identical to the methodology used for the earlier EMEP oxidant model (Jonson et al., 2001). J-values are calculated for clear sky conditions and for two predefined clouds using the PHODIS routine (Kylling et al., 1998). Ozone concentrations from a 2-D global model, extending from the surface to 50 km (Stordal et al., 1985) are scaled by observed total ozone columns from Dutsch (1974). Cloud base for both the predefined clouds is at 1 km above the ground. The first predefined cloud is 3 km deep, with a water content of 0.7 g cm −3 and a mean droplet radius of 10 µm. The second predefined cloud is 1 km deep, with water content of only 0.3 g cm −3 and a mean droplet radius of 10 µm. The J-values are calculated using the recommendations for absorption cross sections and quantum yields from DeMore et al. (1997). The introduction of different chemical mechanisms into the model with new species and photochemical reactions would, in principle, require the recalculation of these databases. As a temporary approach (prior to recalculation of the photolysis databases), we selected the existing photolysis process in the photolysis database which most closely matched the zenith angle dependence of the "new" photolysis process and derived factors to scale the rates. For example, the photolysis of NO 2 provided an excellent description of the photolysis rate of the newly added species, HONO. This is described further in Hayman et al. (2012).

Sulphate production
The parameterization outlined below is previously described in Jonson et al. (2000) and Simpson et al. (2003a), with the only major change being the introduction of explicit pH calculations. In the model SO 2 is oxidized to sulphate both in the gas phase and in the aqueous phase. We always assume equilibrium between gas and aqueous phase. It should be noted that in case the clouds occupy only a fraction of the grid volume, the total concentration (gas + aqueous) of soluble components are assumed to be uniformly distributed in the grid volume. If the cloud evaporates, ions formed in the cloud (e.g. sulphate) are simply assigned to the airborne phase.
For both gas and aqueous phase reactions we scale the reaction rates, rather than the concentrations, by the solubility and cloud volume fractions. In the present calculations we have assumed a constant value cloud liquid water content of 0.6 g m −3 (inside the clouds).
As of model version rv3.9, [H + ] and pH in cloud water is estimated from the acid-base balance, including buffering by bicarbonate (through CO 2 ): This calculation is performed iteratively because the solubility or/and dissociation of SO 2 and NH 3 (and CO 2 ) depend on pH. (Prior to this version a constant pH of 4.3 was assumed). The effect of sea-salt and dust on the cloud pH is not taken into account. Although this could easily be implemented in the model code, large uncertainties are associated with especially the calculations of dust. In any case, studies over continents (and especially industrial/agricultural areas)   Table 5; b for semi-volatile compounds associated with organic aerosol (OA), these characteristics are applied to the particle fraction only.
show that over land cloud-water is dominated by sulphate and nitrate ions and ammonium and hydrogen cations (Aleksic et al., 2009;Aneja and Kim, 1993;Li and Aneja, 1992). The results suggested that the cloud water acidity may be coming predominantly from sulphate aerosol and less from nitric acid. Therefore we have chosen to omit sea salt and dust from the pH calculations. Nitrate and sulphate aerosols and HNO 3 are assumed to be completely dissolved. In the parameterization of aqueous phase chemistry we assume that Henry's law is fulfilled: where [C (aq) ] is the concentration of any soluble gas C (mol l −1 ) in the aqueous phase, H c its Henry's law coefficient and P c the partial pressure of C in the gas phase. In the aqueous phase many soluble gases undergo rapid reversible reactions such as acid-base equilibrium reactions. For these gases it is convenient to define an efficient Henry's law coefficient where the total amount of dissolved gases is taken into ac-count. For example, the total amount of dissolved sulphur in solution (S(IV)) is equal to The total dissolved S(IV) can be related to the partial pressure of SO 2 over the solution (P SO 2 ) by where H SO 2 is the Henry's law coefficient for SO 2 and K 1 and K 2 are the first and second ionisation constants for sulphurous acid.
We define the effective Henry's law coefficient for SO 2 as: and make use of the ideal gas law (P c = [C (g) ] · RT , where [C (g) ] is gas phase concentration of C, R is the universal gas constant and T is temperature) in order to find an expression for the total concentration [C T ] (gas + aqueous-phase) in a cloud volume: where α is the volume fraction of cloud water. Both [C T ] and [C (g) ] are in units M (mol l −1 ). The fraction of the total (gas + aqueous) mass remaining in the interstitial cloud air (f g ) and the fraction absorbed by the droplets (f aq ) can be calculated as: In the model we use the local cloud fraction, defined in the meteorological input fields, as an approximate value for the fractional cloud volume. With the parameterisation above, SO 2 is oxidized both in the cloud free parts of the grid box and in the interstitial cloud air.

Gas phase
In the gas phase SO 2 is oxidized by a chain of reactions initiated by the reaction with OH: with a reaction rate of 2.0×10 −12 cm 3 molecule −1 s −1 . Since some of the SO 2 in a grid square is dissolved in clouds, we define a pseudo reaction rate to allow for this. Using f aq as defined above, then for a fractional cloud volume W , the fraction of SO 2 in the gas-phase is given by: The pseudo-rate coefficient for model reaction OH + SO 2 → SO 4 + HO 2 then becomes k cl−OH = 2.0 × 10 −12 F g (Table S11).

Aqueous phase
Although a number of oxidants may contribute in the oxidation, only O 3 , H 2 O 2 and O 2 catalyzed by metal ions are considered here. The rate of production for sulphate in solution is expressed as: where the reaction rate for the oxidation by O 3 is k cl2 = 1.8 × 10 4 [H + ] −0.4 mol −1 l (Möller, 1980) and the reaction rate for the oxidation by H 2 O 2 is k cl1 = 8.3 × 10 5 mol −1 l (Martin and Damschen, 1981). For the oxidation by O 2 catalyzed by metal ions we assume a reaction rate of k cl3 = 3.3 × 10 −10 s −1 .
As for the gas phase production of sulphate, we define pseudo reaction rates, taking into account the solubility of SO 2 , H 2 O 2 and O 3 and the fractional cloud volume. The pseudo reaction rates then becomes: for the for oxidation by H 2 O 2 , O 3 and O 2 , respectively. f H , f SO 2 and f O 3 are the fractional solubilities of H 2 O 2 , SO 2 and O 3 and is a conversion factor converting k cl1 and k cl2 to molecules −1 s −1 cm 3 . H SO 2 is the Henry's law constant for SO 2 and H * SO 2 is the effective Henry's law constant for S(IV).

Nitrate formation
An important source of aerosol nitrate in the troposphere (and also of NO x loss) is the reaction of N 2 O 5 on deliquescent aerosols, producing two HNO 3 molecules: HNO 3 formed in the reaction above is initially assumed to evaporate and will take part in the formation of ammonium nitrate (Sect. 7.6) or coarse nitrate through reaction 'IN-19' (Supplement , Table S11). Mentel et al. (1999) showed that the uptake rate of N 2 O 5 is around one magnitude lower for nitrate aerosols compared to sulphate aerosols, and this was the basis for the parameterisation of Riemer et al. (2003). More recent measurements in both the laboratory and ambient samples have shown very different values, however, with some studies revealing very low rates, and with very different dependencies, for example on the sulphate/organic ratio (e.g. Brown et al., 2009Brown et al., , 2006Riemer et al., 2009;Chang et al., 2011). Tests with updated schemes have so far not improved the performance of the model for particulate nitrate, and this aspect of the chemistry is probably one of the most uncertain. This reaction is applied whenever RH exceeds 40 %, and following Riemer et al. (2003) the rate we then use is: where c N 2 O 5 is the mean molecular speed for N 2 O 5 and S is here the available aerosol surface area, and α N 2 O 5 is the reaction probability, which is weighted according to the composition of the aerosol: are the aerosol mass concentrations of the secondary inorganic aerosols sulphate and nitrate. (Ideally we would use just fine nitrate here, but given the difficulties associated with such partitioning, we use the more robust sum of fine+coarse nitrate.) The aerosol surface area, S, is calculated from secondary inorganic aerosol mass, m SIA = m SO 2− 4 +m NO − 3 +m NH 2+ 4 , assuming an aerosol density of ρ aer to get volume V , then assuming a log normal size distribution, we get (e.g. Seinfeld and Pandis, 1998): where r n g is the geometric number mean radius (assumed to be 0.068 µm), and σ g = 1.8. The above formulations ignore two terms: (i) the effects of OM and other fine PM on aerosol surface area, which would increase the surface area and hence the rate (ii) inhibiting effect of OM on the sticking coefficient, which would reduce the rate (Riemer et al., 2009). Both terms are very uncertain, but opposite in sign.
For ρ aer we assume a specific aerosol density of 2 g cm −3 near 40 % RH, appropriate for dry aerosol. At higher relative humidity, the salts undergo deliquescence, water content increases, and the density decreases towards values near 1 g cm −3 . The particles grow by absorbing water and hence the surface available to heterogeneous reactions increases. To account in a simple way for the increased surface area, we apply ρ aer = 2.5−1.25RH 100 , RH > 40 (48) where RH is given in %.

Gas/aerosol partitioning
As of version rv3.9, the EMEP model uses the MARS equilibrium module of Binkowski and Shankar (1995) to calculate the partitioning between gas and fine-mode aerosol phase in the system of SO 2+ 4 -HNO 3 -NO − 3 -NH 3 -NH + 4 . MARS has now replaced another code, EQSAM (Metzger et al., 2002;Metzger, 2000), which we have used previously. The MARS module also calculates the mass of aerosol water, see Sect. 11.4.
It should be noted that MARS does not treat sodium chloride and dust components, which is a weakness where seasalt (near coasts) and dust are important. Further, calculated PM water is expected to be underpredicted over seas and coast areas, where sea salt contributes considerably to PM. The effect of not accounting for mineral components is, however, anticipated to be smaller due to their smaller solubility compared to sea salt. It should be recognised that there are also significant uncertainties with other approaches, but in future we will likely replace MARS with a more comprehensive module. (See also Sect. 12.)

Organic aerosol, SOA modelling
As of 2011, a so-called volatility basis set (VBS) approach (Robinson et al., 2007;Donahue et al., 2009) for secondary organic aerosol (SOA) has been added to the available defaults of the EMEP chemical code. The new EmChem09soa scheme uses a variant of the VBS approach which is a somewhat simplified version of the mechanisms discussed in detail in Bergström et al. (2012).
The main differences to the schemes in Bergström et al. (2012) is that in EmChem09soa all primary organic aerosol (POA) emissions are treated as nonvolatile, to keep emission totals of both PM and VOC components the same as in the official emission inventories, while the semi-volatile ASOA and BSOA species are assumed to oxidise (age) in the atmosphere by OH-reactions, leading to decreased volatilities for the SOA, and thereby increased partitioning to the particle phase. We denote this version of the EMEP VBS schemes the "NPAS" scheme (no partitioning of POA, aging of SOA). This assumption of non-volatility for POA is a simplification, but we believe a valid one for our purposes. This is discussed further in Sect. 12.
The OH-reaction rate for SOA-aging in this NPAS scheme is set to 4.0×10 −12 cm 3 molecule −1 s −1 (as suggested by Lane et al. 2008) and each reaction leads to an order of magnitude decrease in volatility and a small increase in mass (+7.5 %) to account for oxygen-addition. This procedure is similar to that used for other EMEP VBS schemes; for further details see Bergström et al. (2012).

Numerical solution of chemical scheme
The chemical equations are solved using the TWOSTEP algorithm tested by Verwer et al. (1996) and Verwer and Simpson (1995). Technical details are discussed in the Supplement, Sect. S2.3.

Resistance formulation
The dry deposition flux (F i g ) of a gas i to the ground surface is modelled using the so-called deposition velocity, V i g (z), such that: This equation is assumed to be true throughout the socalled constant flux layer. In the model we assume that the concentration and deposition velocity calculated at the centre of the lowest grid cell (typically 45 m), a height we refer to below as the reference height z ref , is within this layer. V i g (z) is calculated using a resistance approach: where R a is the aerodynamic resistance between the height z and the top of the vegetation canopy (formally, d + z 0 , where d is the displacement height and z 0 the roughness length), R i b is the quasi-laminar layer resistance to gas i, and R i c is the surface (canopy) resistance.
Over grid-cells which are 100 % sea we simply use the NWP model's meteorological parameters (and z 0 ) to calculate the resistances of Eq. (50). Where grid-cells contain other land-classes, we implement a so-called mosaic approach, whereby the the grid-average deposition rate is given by: where Q symbolises the grid-square average of any quantity Q, f k is the fraction of land-cover type k in the grid-square, and V i g,k is the deposition velocity for this land-cover type, calculated with Eq. (50) using sub-grid (mosaic) values for each resistance term.
In order to make this sub-grid estimation, we are implicitly assuming that the height z ref can be treated as a so-called blending height (e.g. Mason, 1988;Claussen, 1995;Salzen et al., 1996), a height at which the concentrations and meteorological variables are representative of the properties of the full grid square, and not of the local underlying landcover. A further assumption is that the effects of the surface roughness layer can be ignored. Studies have shown that this approximation is probably fine for most purposes, but may impact the estimates of some metrics (AOT40, POD Y , see Sect. 11) (Tuovinen and Simpson, 2008).

Aerodynamic resistance, R a
The first steps in the derivation of sub-grid R a are to derive a grid-square average Obukhov length, L, as in Eq. (8).
The 3-D model meteorology includes wind-speed V H (z ref ) for the centre of the lowest grid level, at around 45 m. We assume that this height is within or near the top of the surface layer, and proceed to calculate turbulence parameters based upon the local values of z 0 and d. These are simply derived from the height, h, of the vegetation for each land-cover type (Table 3). For forests we use d = 0.78h, z 0 = 0.07h, following Jarvis et al. (1976), but with the restriction that z 0 ≤ 1 m. This restriction was found necessary when comparing modelled friction velocity (u * ) with data from the Carbo-Europe network (Papale et al., 2006). For other vegetation, we use d = 0.7h, z 0 = 0.1h. Over water, we use the Charnock relation with z 0 = mu 2 * /g, setting the constant m to be 0.0144 (Garratt, 1992). A minimum value of z 0 = 1.5 × 10 −5 m is enforced, following Berge (1990). From the local d and z 0 , we then estimate a new u * based upon our reference height wind, V H (z ref ): where m is the standard integral function of the similarity profile of momentum (Garratt, 1992). Having calculated u * in this way, a local estimate of L can be found by substituting u * in Eq. (8). The aerodynamic resistance for heat or scalars between any two levels z 1 and z 2 is calculated with the standard R a (z) formula, the same as used in Eq. (50).

Quasi-laminar layer resistance, R i b
The quasi-laminar layer resistance is calculated with where Sc i , the Schmidt number is equal to the ν/D i , with ν being the kinetic viscosity of air (0.15 cm 2 s −1 at 20 • C) and D i is the molecular diffusivity of gas i, and P r is the Prandtl number (0.72). Over sea areas the expression of Hicks and Liss (1976) is used:

Surface resistance, R c
Surface (or canopy) resistance is the most complex variable in the deposition model, as it depends heavily on surface characteristics and the chemical characteristics of the depositing gas. Our approach makes use of bulk canopy resistances and conductances (R and G terms, where G i = 1/R i for any gas i), and of unit-leaf-area (one-sided, projected) resistances and conductances, which we denote with lowercase letters (r, g). The general formula for bulk canopy conductances, G c , is: where LAI is the one-sided (projected) leaf-area index (m 2 m −2 ), g sto is the stomatal conductance, and G ns is the bulk non-stomatal conductance. For non-vegetative surfaces only the last term is relevant. The formulation for stomatal and non-stomatal conductances for most gases and conditions are dealt with in Sects. 8.5-8.6. Two special cases are (a) HNO 3 and (b) NH 3 over crops:

(a) R c , HNO 3
In normal conditions the surface resistance to HNO 3 is effectively zero. A minimum value of R c of 10 s m −1 is enforced for numerical reasons, so for HNO 3 the whole canopy resistance is then simply given by: where R HNO 3 low accounts for observations of HNO 3 deposition over snow, and is set simply to R HNO 3 low = −2 T S , with T s being the surface (2 m) temperature in • C. These values loosely match those found by Johansson and Granat (1986) for temperatures of down to −18 • C. During the growing season for crop land-covers, the surface resistance is set very large, ensuring zero deposition. This procedure is designed to account for the fact that many croplands are actually emitters of NH 3 , rather than sinks (e.g. Sutton et al., 2000;Fowler et al., 2009, and references therein).

Stomatal conductance, g sto
Stomatal conductance is calculated with a multiplicative model, a development of that described in Emberson et al. (2000a): where g max is the maximum stomatal conductance, and f x are factors (within 0-1) accounting for time of year (leaf phenology), the minimum observed stomatal conductance (min), light (actually PAR), temperature (T ), vapourpressure deficit (D), and soil-water (SW). It should be noted that the canopy scale stomatal conductance (LAI g sto in Eq. 55) is a non-linear function of LAI, since f light and hence g sto are non-linear functions of LAI, see Supplement, Sect. S7.2. The main new feature of the EMEP model with regard to this procedure is that soil water effects are now included by default. In Emberson et al. (2000a), f SW was based upon soil-water-potential (SWP). SWP is a very non-linear function of soil water content, varying with soil texture and homogeneity, and in practice can only be accurately estimated with in situ measurements. For these reasons f SW was simply set to 1 in most previous EMEP model runs, i.e. stomatal uptake was not assumed to be limited by soil water availability (e.g. Simpson et al., 2007a). A number of techniques are being investigated with regard to soil water calculations (Büker et al., 2011), but as of version rv3.9 the EMEP code makes use of the simple S MI index (Sect. 3.3) to calculate f SW . Rather than attempting to calculate absolute values of SWP, we use a simple procedure designed to capture the main effects of dry periods on g sto : For deposition modelling we use the S MI values appropriate for deeper soil layers; for ECMWF inputs this is the top 1 m soil layer. The methodology for g sto was developed and tested within a dry deposition framework for ozone, now referred to as the DO 3 SE (Deposition of Ozone and Stomatal Exchange) model (Emberson et al., 2000a(Emberson et al., ,b, 2001(Emberson et al., , 2007Klingberg et al., 2008;Simpson et al., 2001Simpson et al., , 2003bTuovinen et al., 2001Tuovinen et al., , 2004. Stomatal conductance calculated for any other gas i is simply scaled from that for ozone using the ratio of the diffusivities in air of ozone and gas i. Table S18 in the Supplement gives the diffusivities (although expressed relative to water) used in the EMEP model.
Further details of the equations and current parameter values underlying the stomatal conductance algorithm are given in the Supplement, Sect. S7.2.

Non-stomatal resistances
G ns is calculated specifically for O 3 , SO 2 , and NH 3 . Values for other gases are obtained by interpolation of the O 3 and SO 2 values (Sect. 8.8).
The ground-surface resistance, R i gs , for a specific gas is an important component of the total non-stomatal resistance.
Base-values of R gs (denotedR gs ) for O 3 or SO 2 are given in Table S19. Similar to Zhang et al. (2003), these are modified for low temperature and snow cover with: where x represents either O 3 or SO 2 , f snow reflects the snow coverage, and F T is a low-temperature correction factor -see Sect. 8.7.1 for both terms.

Ozone, G O 3 ns
Our formulation of the non-stomatal conductance for ozone builds upon the framework of Emberson et al. (2000a), which has been been extensively evaluated in a number of studies (Emberson et al., 2000a;Tuovinen et al., 2001Tuovinen et al., , 2004: where SAI is a surface area index (m 2 m −2 ), r ext is the external leaf-resistance (cuticles+other surfaces) per m 2 PLA, R inc is the in-canopy resistance, and R gs is the ground surface resistance (soil or other ground cover, e.g. moss). The external resistance r ext is set to 2500 F T s m −1 , where F T is a low-temperature correction factor (see Sect. 8.7.1). Following Erisman et al. (1994), the in-canopy resistance, R inc , is defined as b SAI h/u * , where h is the canopy height and b = 14 s −1 is an empirical constant. SAI is simply set to LAI+1 for forests, or equal to LAI for non-crop vegetation. For crops a substantial part of the leaf area can be senescent. A simplified version of the methodology of Tuovinen et al. (2004), based upon the life-cycle of wheat, is applied: where d N is the day number, and d SGS , d EGS , and L S are as defined in Sect. 5. Outside the growing season, SAI = LAI = 0. For vegetated surfaces, the non-stomatal resistance R ns for NH 3 is assumed to depend upon surface (2 m) temperature, T s ( • C), humidity levels, RH (%), and on the molar "acidity ratio":

Atmos
This acidity ratio is a first attempt to account for the observed changes in resistance in areas with different pollution climates Fowler and Erisman, 2003). More advanced treatments are possible, but the spread in values from different parameterisations is substantial (Massad et al., 2010).
The parameterisation of Smith et al. (2000) has been modified in order to take into account the effects of a SN , based upon an approach suggested by Smith et al. (2003). The resulting scheme can be expressed as: where β is a normalising factor (1/22 = 0.0455). The F 1 term is identical to that of Smith et al. (2000) and provides a relationship of R ns with temperature and relative humidity. The second function, F 2 , is an equation derived from observations presented in , and relates the value at 95 % relative humidity and 10 • C to the molar ratio of SO 2 /NH 3 . The two terms are equal for molar SO 2 /NH 3 ratio 0.3. The factor β is introduced in order to normalize one equation to the other, i.e. to ensure that the combined parameterisation is equal to the two separate terms for 95 % relative humidity, 10 • C and molar ratio 0.3. For above-zero temperatures R NH 3 ns is constrained to lie between 10 and 200 s m −1 . Finally, we do not distinguish wet or dry surfaces in this formulation (they are included in the RH dependency used above).

SO 2 ns
The canopy conductance of SO 2 is strongly controlled by wetness and NH 3 levels, as well as deposition of other acidic gases (HNO 3 and HCl), adsorption of CO 2 , aerosol dry deposition, the composition of rain during precipitation events, ion leaching from the plants and processes such as dew fall and guttation (e.g. Flechard et al., 1999;Fowler et al., 2001Fowler et al., , 2009Burkhardt et al., 2009).
In order to develop a simple parametrisation for G SO 2 ns , which nevertheless captured the main processes, Fagerli et al. (2012) used long-term simultaneous measurements of NH 3 and SO 2 exchange, made within the EU LIFE Deposition Monitoring Project , to derive operational parameterisations of co-deposition effects.
The parameterisation developed links the non-stomatal canopy uptake resistance of SO 2 (R SO 2 ns ) to the mean molar SO 2 /NH 3 ratio in air over the last 24 h, a 24h SN : For above-zero temperatures R SO 2 ns is constrained to lie between 10 and 1000 s m −1 . a 24h SN is constrained to be maximum 3, which corresponds to R SO 2 ns = 400 s m −1 for RH of about 85 %. For non-vegetative surfaces, R SO 2 ns is simply set to the base-values,R gs , shown in the Supplement, Table S19.

Snow and low-temperature corrections
At temperatures below −1 • C, non-stomatal resistances are increased using a factor F T as in Zhang et al. (2003): with the constraint 1 ≤ F T ≤ 2. Resistances for SO 2 over snow covered surfaces depend on the temperature. For instance, Granat and Johansson (1983) found that SO 2 dry deposition velocities were smaller than 0.1 cm s −1 at temperatures below −1 • C, but higher at warmer temperatures due to the presence of liquid water at the snow surface. R snow for SO 2 (in s m −1 ) are here loosely based on Erisman et al. (1994) and Zhang et al. (2003): For ozone, we simply set R O 3 snow = 2000 s m −1 . The term f snow in Eq. (59) is an estimate of the fractional cover of snow, derived from the NWP model's snow depth (S d ) and an assumed maximum value S d,max at which the snow fraction for canopy leaves is assumed to be 1. We use a similar methodology to that proposed by Zhang et al. (2003): with the constraint 0 ≤ f snow ≤ 1. Zhang et al. (2003) presented tabulated values of S d,max , but we simply assume that S d,max = 0.1 h, where h is the height of the vegetation. If some fraction of the grid is covered with ice, we assume that f snow is the maximum value of the snow or ice fractions.

Extension to other gases
For all gases other than HNO 3 or NH 3 we obtain G ns by interpolating between the values for O 3 and SO 2 . This interpolation borrows the solubility index, here denoted H * , and the reactivity index, f 0 , from the Wesely (1989) methodology, but these are applied directly now to total non-stomatal conductance rather than to individual resistances (Table S18). As there is so little data available on non-stomatal resistances, even for O 3 and SO 2 , this simpler scaling seems acceptable. With these indices, the dry and wet conductance values for a gas i are obtained from the values for ozone and SO 2 using:

Aerosol dry deposition
Although a range of theory-based models is available to describe aerosol deposition, they often predict features which conflict with measured deposition rates (Pryor et al., 2008b,a;Petroff et al., 2008a;Flechard et al., 2011). For example, methods based on the well-known formulations of Slinn (1982) predict low deposition velocities to forest canopies. Alternative formulae of Zhang et al. (2001) predict higher deposition velocities, but no effect of canopy density. Several studies show that ammonium-nitrate has higher deposition velocities than sulphates, as a result of the partitioning of NH 4 NO 3 to the more rapidly depositing HNO 3 and NH 3 gases (e.g. Fowler et al., 2009;Wolff et al., 2010). Petroff et al. (2008a,b) have presented an extensive discussion of the issues surrounding chemically-inert particles, and presented calculations where deposition is affected by both particle size and canopy leaf area index. Loosely based upon these reviews, and results from various experimental studies, we have implemented a new but deliberately simple scheme for particles in low vegetation and forests in the EMEP model. The basic formulation follows the same pattern as many studies (Wesely et al., 1985;Lamaud et al., 1994;Gallagher et al., 1997;Nemitz et al., 2004), but modified by an enhancement factor, F N , for nitrogen compounds: where V ds is the surface deposition velocity (Petroff et al., 2008a), and F N = 3 for fine-nitrate and ammonium, and 1 for all other compounds (Table 6). Further, we restrict application of the equation to 1/L > −0.04 m −1 . For all landcover categories except forests we use use a 1 = 0.002 from Wesely et al. (1985), and set a 2 to 300 m, the simplified stability correction suggested by Gallagher et al. (1997).
For forests, we implement a simple dependence on surface area index: with a 2 again set to 300 m, and the additional restriction that a 1 ≥ 0.002. These values are loosely based upon the results of an analysis of measurements, and sets of complex calculations presented in Petroff et al. (2008a,b). Petroff et al. (2008b) calculated that a forest with total LAI of 22 would have a surface deposition velocity of ca. 0.002-0.004 m s −1 at u * = 0.45 m s −1 for particles in the accumulation size range (see Fig. 15, Petroff et al. 2008b). Our 0.008 u * gives 0.004 m s −1 for this same friction velocity. They also showed that a decrease in LAI of a factor of 2 would reduce V ds by a factor 1.5-2. Further, Petroff et al. (2008b)'s calculations suggested that V ds is approximately proportional to LAI for D p ∼ 0.5 µm. For the EMEP model we make use of our surface area index, SAI, which accounts for non-leafy surfaces, and which is simply derived as SAI = LAI + 1 for forests. Petroff started with a total LAI of 22, which is ca. LAI = 10 (1-sided), or SAI = 11. Simplifying, we therefore scale with SAI/10. (The use of SAI rather than LAI also prevent wintertime deposition in deciduous forests going to zero). Finally, we enforce a minimum V ds of 0.002 u * , consistent with Wesely as SAI → 0.
As pointed out by Venkatram and Pleim (1999), the resistance analogy is not appropriate for particles. We have therefore implemented the mass-conservative equation: where v s is settling velocity, V d (z) is the deposition velocity at height z, and r(z) is the sum of the aerodynamic resistance and inverse V ds . As summarized in Sect. 6, the EMEP model distinguishes five classes of fine and coarse particles, which are presently assigned mass-median diameters, geometric standard deviations (σ g ), and densities (ρ p ).
Although the dry-deposition rates of fine (accumulationmode) particles are not size-dependent in the model, the overall dry deposition rate of larger particles is affected by v s , which is strongly size-dependent. To account for this, the v s calculations are integrated over the aerosol sizes, assuming a log-normal particle size distribution. These polydisperse settling velocities of coarse particles are calculated, using Eqs. (A25-A32) from Binkowski and Shankar (1995).
This revised scheme (and the changes in assumed aerosol size), which we here denote the EMEP-12 particle deposition scheme, gives significantly different rates to those used previously, with higher rates for fine particles (especially for the nitrogen components), and lower rates for coarse nitrate (since the assumed particle size is smaller). In order to illustrate the net effect, and place these results in the context of a previous comparison, we have rerun the setup of Flechard et al. (2011), but adding the new EMEP-12 particle deposition scheme. In Flechard et al. (2011), four different deposition modules (including EMEP-03) were applied for 55 European sites covering four land-cover categories: Forest (F), Seminatural (SNL), grassland (G) and crops (C). This study also made the general assumption that 19% of nitrate is in the coarse mode at all sites. The sites were part of the EU Ni-troEurope study, monitoring monthly concentrations of the key reactive nitrogen (Nr) species, with the intention to estimate dry-deposition using inferential techniques (Tang et al., 2009). Figure 6 compares the estimated deposition rates for particulate nitrate from the 2003 and 2012 versions of the EMEP scheme, and three other models, CBED, CDRY and IDEM models (for references and descriptions of these other models, see Flechard et al. 2011). The EMEP-03 and EMEP-12 results are surprisingly similar for all land-cover categories except crops, where EMEP-12 gives higher rates. This similarity is partly coincidental, however, representing the balance between increased deposition rates for fine particles, and reduced rates for coarse particles. For example, over grassland the estimated V g for fine-nitrate increased by a factor of 5 on average (from 0.28 mm s −1 for EMEP-03 to 1.4 mm s −1 for EMEP-12), but V g for coarse nitrate decreased by a factor of six (from 7.2 to 1.2 mm s −1 ). The larger change for crops seen in Fig. 6 reflects the more complicated changes in R b (with different equations used inside and outside the growing season) used in the EMEP-03 scheme. The changes from EMEP-03 to EMEP-12 are thus significant, but as also seen in Fig. 6, differences between all methods are very large. As noted in Flechard et al. (2011), this is unfortunate, but currently the experimental difficulties are too large to allow a reliable choice of scheme (e.g Fowler et al., 2009;Pryor et al., 2008b,a). The new EMEP particle deposition scheme has at least the advantages of simplicity of formulation, and results are broadly consistent with recent but more complex schemes, and recent flux measurements Petroff et al., 2008a,b;Wolff et al., 2010).

Wet deposition
Parameterisation of the wet deposition processes in the EMEP model includes both in-cloud and sub-cloud scavenging of gases and particles. The parameterization of the wet deposition is previously described in Berge and Jakobsen (1998).

In-cloud scavenging
The in-cloud scavenging S in of a soluble component of mixing ratio χ is given by the expression: where W in is the in-cloud scavenging ratio given in the Supplement, Table S20, P (kg m −2 s −1 ) is the precipitation rate, h s is the characteristic scavenging depth (assumed to be 1000 m) and ρ w is the water density (1000 kg m −3 ). We do not account for the effect that dissolved material may be released if clouds or rain water evaporate.

Below-cloud scavenging
For below cloud scavenging a distinction is made between scavenging of particulate matter and gas phase components. The sub-cloud scavenging of the gases is calculated as: where W sub is the sub-cloud scavenging ratio given in the Supplement, Table S20. Wet deposition rates for particles are calculated, based on Scott (1979), as: where V dr is the the raindrop fall speed (V dr = 5 m s −1 ), A = 5.2 m 3 kg −1 s −1 is the empirical coefficient (a Marshall-Palmer size distribution is assumed for rain drops), andĒ is the size-dependent collection efficiency of aerosols by the raindrops (Table S20). The collection efficiency is size dependent, with a minimum for fine particles (see Laakso et al., 2003;Henzing et al., 2006).
Initial concentrations of major long-lived species are required in order to initialise model runs. Boundary conditions along the sides of the model domain and at the top of the domain are then required as the model is running. Additionally, we often need to specify concentrations of some species which are not explicitly included in the chemistry of interest, but that enter into reactions with some of the reacting chemical compounds ("background" species). We refer here to all of these types of data as initial and boundary conditions (IBCs). Two main methods of specifying boundary conditions are currently available: 1. Provision of 3-D fields for whole domain from previous runs of the same or another version of the EMEP model (self-assimilation), or from other models, typically global chemical transport models (CTMs).
2. Simple functions are used to prescribe concentrations in terms of latitude and time-of-year, or time-of-day. For ozone, 3-D fields for the whole domain are specified from climatological ozone-sonde data-sets, modified monthly against clean-air surface observations.
Method (1) allows great flexibility. A pre-processing program interpolates the data field of interest to the desired horizontal resolution (e.g. 50 km × 50 km), and to the 20 vertical levels in the EMEP model. The frequency of the update of the boundary conditions can be chosen freely, as long as the boundary condition field is provided for the same time period. Examples of this kind of approach can be found in Vieno et al. (2010), where the European scale model was used to provide IBCs for a 5 km scale model over the UK.
Method (2) is used for those species where rather simple descriptions of boundary condition are sufficient. Despite its simplicity, this method has the advantage that the IBCs are based upon measurements, ensuring a robustness which global CTM model results sometimes lack. For policy runs, the EMEP model is usually run using this methodology, and it is this method we document here.

Ozone
Ozone is the gas for which the specification of accurate boundary conditions is most essential to a good model performance. This is due to the fact that ambient ozone levels in Europe are typically not much greater than the Northern hemispheric background ozone. Boundary conditions of ozone are developed from a two-step procedure. First, the climatological O 3 data of Logan (1998) is used, which provides gridded O 3 data with resolution 4 • latitude by 5 • longitude for 13 pressure levels. These data are interpolated to the EMEP grid system to provide a monthly base-set for ozone IBCs.
These monthly data are then adjusted using a so-called "Mace-Head" adjustment. Mace Head is a site on the west coast of Ireland, ideally suited as a background site for midlatitude air masses. It was shown by Derwent et al. (1998), using trajectory analysis and other techniques, that the cleanair concentrations of O 3 (and CO) at Mace Head were basically uniform in a wide sector for air masses arriving from Iceland to Barbados -in other words, it confirmed the view of a general well-mixed background air mass.
For the EMEP model we have made use of an extended version of this analysis. Ozone concentrations from Mace Head have been sorted using sector-analysis, obtained using trajectories obtained from http://www.emep.int 1 . Monthly mean values of the ozone associated with easterly sectors (sectors 6-8) have been calculated. Where fewer than 15 days were available to make an average for a particular year, averages from a full 10-yr analysis were substituted for the missing days.
In order to generate an adjustment factor, the monthly values of observed O 3 derived using this procedure, is the average concentration from model domain x = 1..x M , y = 1..y M ). If the difference between the two datasets obtained in this way is (=O MH 3 -O GD 3 , in ppb), we simply add to the ozone boundary conditions over the whole domain. Since the concentrations of ozone are generally increasing with height in the model domain (from say 40 ppb to several hundred ppb), then the effect of this constant term is greatest for the surface layer and quite small at say 5-10 km height.
Although simple, this procedure ensures that the BCs used for ozone are realistic in the mid-latitude region near ground level, at least near the Western boundary. Although based entirely upon one station, this correction has been found to result in good BCs for almost all sites on the west coast of Europe, ranging from Norway to Spain.
For other species where prescribed values are needed, simple functions have been chosen, designed to enable concentration values that correspond to observations. The concentrations are adjusted in the vertical and for latitude and time of the year (monthly fields) to match the observed distributions. Table S21 lists the parameters used, as described below.
We first calculate the seasonal changes in ground-level BC concentration, χ 0 , through: where χ mean is the annual mean near-surface concentration, χ the amplitude of the cycle, n y is the number of days per year, d mm is the day number of mid-month (assumed to be the 15th), and d max is day number at which χ 0 maximises. Changes in the vertical are specified with a scale-height, H z : where χ IBC (h) is the concentration used for IBCs at height z.
For simplicity we set z to be the height of the centre of each model layer assuming a standard atmosphere. Values of χ IBC are constrained to be greater or equal to the minimum values, χ v min , given in Table S21. For some species a latitude factor, given in Table S22, is also applied. Values of χ i adjusted in this manner are constrained to be greater or equal to the minimum values, χ h min , given in Table S21. Finally, for two species, we simply specify constant mixing ratios over the whole model domain, valid for 1990 (see Sect. 10.2 for other years). These are 1780 ppb for methane and 600 ppb for hydrogen.

Trends in initial and boundary conditions
The BC values discussed above are assumed appropriate for the year 1990. For other years these values are adjusted using trend factors. Such adjustments can be made with results of e.g. global CTMs (including EMEP model runs). Lacking other information we use the default trend factors as summarised in the Supplement, Table S23.

Outputs
The EMEP model produces a large number of outputs for a variety of purposes. Most are straightforward, for example maps of annual wet deposition of oxidised or reduced nitrogen. However, some outputs display special features or are provided for specific purposes. For example, one of the main reasons for running the EMEP model is to generate results for use in integrated assessment modelling (IAM), and for studies on the risks and damages caused by pollution, and a number of model outputs are designed with this in mind.
Here we briefly describe some of the most important outputs.

Near-surface concentrations
The basic calculations of the EMEP CTM produce concentrations for model layers. The lowest layer is about 90 m deep, so concentrations from this layer may be interpreted as being applicable for 45 m above ground level (or stricter, above displacement height d). In order to estimate concentrations at heights more typical of measurements, typically around 3 m for EMEP observations, or at canopy top for some ozone-flux or AOT40 estimates, we make use of assumption that the vertical deposition flux density (F i g , Eq. 49) remains approximately constant within the atmospheric surface layer (e.g. Tuovinen, 2000). Referring to the model concentrations of species i at reference height z ref of 45 m as χ i (z ref ), we readily obtain the concentrations at any other height within the surface layer from Eq. (49): with appropriate calculations of the deposition velocity resistance terms as discussed in Sect. 8.

Ecosystem-specific depositions
As discussed in Sect. 8, the model's calculations of dry deposition are made separately for each sub-grid landcover. For provision to IAM or the effects community, these sub-grid estimates are aggregated to provide output deposition estimates for broader ecosystem categories, as shown in Table 7 A possible output would be deposition to water, but for IAM purposes the deposition of interest here is to the catchment area, rather than to the water surface. Thus, deposition estimates for waters are usually simply taken from the gridaverage depositions.

Ozone statistics
A number of statistics are typically used to describe the distribution of ozone within each grid square, and for input to IAM assessments: Mean of Daily Max. Ozone. -First we evaluate the maximum modelled concentration for each day, then we take either 6-monthly (1 April-30 September) or annual averages of these values.

AOT40 uc
c . -AOT40 calculated for agricultural crops using estimates of O 3 at the top of the crop. This AOT40 is close to that defined for agricultural crops by LRTAP (2009), but using a default growing season of May-July, and a default crop-height of 1 m.
-as above, but using the simple grid-average concentrations from the model's 3 m level.
The first two "canopy-top" definitions are in accordance with the recommendations of LRTAP (2009), and the two "grid" values are for comparison to AOT40 maps derived from observations. In all cases only daylight hours are included, and for practical reasons we define daylight for the model outputs as the time when the solar zenith angle is equal to or less than 89 • . (The proper UNECE definition uses clear-sky global radiation exceeding 50 W m −2 to define daylight). The EU definitions of AOT40 use day hours from 08:00-20:00.
For the development of the 1999 "Gothenburg" Protocol (http://www.unece.org/env/lrtap/), the metric used for assessing the risk to vegetation was AOT40. However, new critical levels based on POD Y have now been agreed (Mills et al., 2011b, and references therein). For provision of data to support the use of these new approaches to IAM, a simplified approach to mapping ozone fluxes was defined by LRTAP (2009), in which one generic crop species was defined, and two generic forest species. The "IAM" species in Tables 3  and Table S16 correspond to these, although the phenology functions are somewhat simplified compared to the latest (2010) Mapping Manual update. In the model inputs, a tiny fraction of IAM CR, IAM DF and IAM MF are added to each grid square where any vegetation is present, so we can calculate fluxes even in grids where the landuse data suggest no such species are present, providing a more comprehensive and easier to interpret spatial indication of risk.
This simplified approach for IAM was adopted because it was recognised that our knowledge of many critical inputs Table 7. Ecosystems provided in deposition outputs, and associated EMEP landcover categories (see Table 3). (e.g. growing seasons and phenology, conductance parameters, elevation effects, soil water parameters, etc.) is too uncertain to allow accurate mapping of the real ozone flux to specific species. On the other hand the spatial distribution of fluxes is so different to that of AOT40 (Simpson et al., 2007a) that calculation of fluxes to a generic species was seen as an improvement upon the continued use of AOT40. It was also recognised that the IAM process (which balances health and vegetation impacts from many pollutants, against costs of emissions measures) could not take into account many different types of vegetation, and that only a few flux-maps could be included in the IAM optimisation work.
Although there are obvious similarities in the methods used to model upper-canopy stomatal fluxes (F st ) for the calculation of POD Y , and modelling of full-canopy fluxes for deposition purposes, these calculations have important differences. The F st values required for POD Y represent maximum uptake to a small portion of the canopy, not net uptake to the whole canopy. These F st calculations are therefore performed as a parallel exercise to the deposition modelling, being performed from within the EMEP model's deposition routines, but having no feedback to the canopyscale deposition calculations required for the model's atmospheric chemistry calculations. The f light term (see Supplement, Sect. A6.2) is based upon I sun PAR , and soil-water limitations usually ignored (i.e. f SW = 1). Further discussion of these type of calculations is given in Simpson et al. (2007a) and Tuovinen et al. (2009).
For these generic "IAM" species, the suffix gen can be applied, e.g. POD Y,gen is used for forests. (POD was introduced in 2009 as an easier and more descriptive term for the accumulated ozone flux than the former AFst term. The definitions of AFst and POD are identical however.)

PM-water
PM 10 and PM 2.5 mass determined with a gravimetric method is likely to include particle-bound water, which does not get completely removed (or condenses on the particles) under filters conditioning at temperature 20 • C and relative humidity 50 %. To make comparison of calculated PM 10 and PM 2.5 concentrations with gravimetric measurements more consistent, the model accounts for particle water within the PM mass. The water content in PM 2.5 and PM 10 is calculated Atmos. Chem. Phys., 12, 7825-7865 with the MARS equilibrium model (Binkowski and Shankar, 1995) for the conditions required for filters equilibration, i.e. temperature 20 • C and relative humidity 50 %. As only fine SIA aerosols (i.e. SO 2− 4 , NO − 3 and NH + 4 ) are included in the MARS model, the calculated water describes water in PM 2.5 . The calculated mass of water is added to both dry PM 2.5 and PM 10 masses when being compared with measured concentrations. Note that the components of sea salt aerosol is not included in the MARS model, leading to some underestimation of particle water.
The calculated aerosol water content depends on the mass of soluble PM fraction and on the type of salt mixture in particles. Accounting for particle water in calculated PM 2.5 and PM 10 has been shown to improve the general correspondence between model results and observations. However, there are caveats to the model esimates of particle-bound water as no proper verification of the calculated water content with measurements is presently available. Further details as well as results and initial evaluation of model calculation of particle water can be found in Tsyro (2005).

Discussion and some future challenges
As noted in the Introduction, the intention of this paper has been to document the EMEP MSC-W model version rv4, and for reasons of space this has not allowed much discussion of the background for model choices. One motivation for this focus on documentation is of course the importance of the EMEP model for European air pollution policy formulations. Another motivation is that the model is being widely used, applied in research projects and/or model intercomparison excercises, but so far only sparse and incomplete documentation has been available for the recent model versions. Indeed, the model has taken part in a large number of inter-comparisons in recent years: Cuvelier et al. (2007); van Loon et al. (2007) 2012), among others. In terms of performance, the EMEP model has ranked well in these studies, with consistently good performance for different pollutants (ozone, PM, etc.). In terms of complexity the EMEP model is fairly similar to other regional-scale European CTMs, such as MATCH (Robertson et al., 1999), CHIMERE (Bessagnet et al., 2004), or DEHM (Christensen, 1997;Frohn et al., 2001). All of these models have some flexibility with regard to chemical schemes, and have zooming-capabilities.
Given the complexity of any CTM, it is hard to limit a discussion of where the main limitations in a model are, and indeed it is difficult to know if the main source of uncertainty in models lies in their meteorological drivers, physical descriptions, chemical and/or aerosol schemes, or loss processes. The reliability of inputs such as emissions is a major cause of uncertainty. Here we address just a few areas where improvements are desired in the next few years, and where some work is ongoing. To limit the scope, we focus on particulate matter, which is probably the biggest challenge for both CTM models and policy development.

(i) Aerosol size-distributions
The standard EMEP model described here uses essentially two size-modes for particles, although our definitions of particle-size depend a little on the compound. This is a great simplification, which can be justified for current needs by the fact that the present version of model is mainly designed to calculate PM 10 and PM 2.5 mass closure (i.e. concentrations and chemical composition), which over the last decade has been the highest priority within the EMEP/LRTAP Convention framework. A pragmatic defence of this procedure is that in most comparisons with measurements (e.g. Fagerli and Aas, 2008;Simpson et al., 2006b, or EMEP status reports over many years) the EMEP model has been shown to perform quite well against measured PM mass. Problems are clearly apparent in some studies, for example in capturing hourly variations in Nr components measured in reactive nitrogen components close to agricultural areas (e.g. Aas et al., 2012), but it is unclear how far this problem can be related to size-distribution, and how much is due to other (so far unsolved) problems with model resolution or equilibrium dynamics (e.g. Aan de Brugh et al., 2012).
Still, the need for a more detailed description of the aerosol size-distribution is clearly apparent, on grounds of scientific realism (to capture the effects of for example in-cloud activation of particular size-fractions), an increasing need to link to climate issues (e.g. Liu et al., 2012), and also in terms of human health effects, where size distributions are also implicated. The challenge here is mainly to find an optimal balance between the number of bins/modes, in order to increase realism but preventing excessive CPU-increases.

(ii) Gas/aerosol equilibria (inorganic)
As noted in Sect. 7.6, the MARS module we use for gas/aerosol partitioning of inorganic compounds into finemode aerosol cannot account for sea-salt and dust components, and we use a very simplified treatment of nitrate formation on coarse aerosol. In future we will likely replace MARS with a more comprehensive module (e.g. Fountoukis and Nenes, 2007), and likely use a kinetic (rather than equilibrium) approach for coarse nitrate formation, with explicit reactions of for example HNO 3 with NaCl or dust. We have indeed been exploring such reactions, but this is ongoing work. Apart from the difficulty of predicting such components, there are also large gaps in our scientific understanding of nitrate composition -there are hardly any measurements of coarse-mode nitrate to compare against for example. For organic aerosol (OA), there are a large number of problems with all model formulations, something which inevitably follows from the complexity of OA itself, and our lack of understanding of the underlying science (e.g. Hallquist et al., 2009;Ng et al., 2010). Many of these issues (e.g. large uncertainties in emission inventories from both anthropogenic and biogenic sources, or in vapour pressure assumptions) have been discussed in relation to earlier EMEP modelling studies (Simpson et al., 2007b;Bergström et al., 2012), but here we briefly discuss an uncertainty arising from the use of the simplified 'NPAS' VBS scheme in the standard model. The NPAS scheme assumes that POA emissions can be treated as non-volatile, instead of treating them (and related emissions of SVOC and IVOC: semi-and intermediatevolatilty gases) as components of varying volatility (as in, e.g Robinson et al., 2007). In high-emission areas this NPAS scheme should lead to higher OA compared to a model that allows evaporation of some of the initially emitted POA. On the other hand, VBS schemes often postulate emissions of SVOC & IVOC which are supposed to be unaccounted for in the official emission inventories (e.g. Shrivastava et al., 2008). This is likely more realistic, and provides a larger pool of VOC compounds from which partitioning to aerosol can occur, but it is also a large source of uncertainty.
There are two main reasons why we choose to use nonvolatile POA emissions in the "standard" EMEP model code (that used for policy-associated runs): (1) The volatility distribution of POA and associated SVOC and IVOC compounds is poorly known; the amount of SVOC + IVOC emissions is probably substantial, but so far we have only a very limited number of (American) studies with which to estimate this contribution (e.g. Shrivastava et al., 2008); (2) official European emission inventories used for policy modelling consist of PM compounds which are assumed to be inert, as well as VOC emissions. No consideration of volatility is made in either the PM or VOC inventories. For policy modelling it is necessary to keep these POA and VOC emission totals the same as in the official emission inventories.
In order to assess the sensitivity of the model to this assumption, we have used the schemes presented by Bergström et al. 2012 to compare model versions with and without this inert POA assumption. The results, illustrated in Supplement Fig. S2, show that total yearly average OA concentrations (in PM 2.5 ) over most of the European land-area are 10-20 % lower when we use inert POA emissions (NPAS scheme) than if we use volatility-based emissions and aging of the emitted semi-and intermediate volatility OC emissions -an effect of the extra SVOC+IVOC assumed in the PAA VBS scheme which generates more OA because of aging processes. For some high-emission areas the inert assumption leads to much higher yearly average OA than the volatility based approach. The biggest effects are found over Paris, where we obtain more than 40 % higher fine OA with the inert POA model.
The volatility question is thus important, but one of many uncertainties with regard to OA modelling. There is clearly an urgent need for new measurement and inventory data on emissions and volatility distributions in Europe though, if we are to take account of these properly in future research or policy-related modelling.

(iv) Dispersion/resolution issues
The final major challenge we will mention is that of model resolution, especially in the vertical. The EMEP model currently has a lowest layer of about 90 m deep, so generates concentration data which is appropriate for about 45 m. This is clearly a simplification, especially in wintertime for those pollutants that have major ground-level sources. Woodburning emissions are a clear example of this, and their importance for wintertime OA concentrations has been stressed in many studies (see Simpson et al., 2007b;Bergström et al., 2012, and references cited therein), but there is no easy solution. Simply reducing the thickness of the lowest model layer is also unrealistic, since this quickly leads to an imbalance in the horizontal and vertical scales. It is not realistic for example to disperse the emissions from a point source through 50 km × 50 km in the horizontal, but only 10 m in the vertical.
Another major problem for this situation is that modelling of the stable boundary layers found in wintertime is notoriously difficult, and there is generally very little data with which to evaluate the vertical dispersion estimated by models. Work is also ongoing to find the best compromise in terms of grid resolution and model physics.

Conclusions
The Meteorological Synthesizing Centre -West (MSC-W) of EMEP has been performing model calculations in support of UNECE for more than 30 years. The EMEP MSC-W chemical transport model is still one of the key models used in policy support in Europe. It is central to UN-ECE work, with a mandate to provide scientific support to the development of air pollution reduction Protocols, and is the sole provider of source-receptor matrices to the IIASA GAINS model (which is central to EU policy work), and is used in many EU projects alongside other chemical transport models.
The MSC-W models have been increasing in complexity and capabilities over this time-period, and today the MSC-W model is used to simulate photo-oxidants and both inorganic and organic aerosols, on scales ranging from national studies at ca. 5 km resolution to global scale.
The last full documentation of the EMEP MSC-W model is almost ten years old (Simpson et al., 2003a), referring to version rv1.7 (or EMEP-03 for simplicity). The model has changed in numerous ways (both large and small) since this document was written. These changes include revised methods for calculating mixing heights and eddy diffusion coefficients (for stable and neutral conditions), new temporal variation factors (based upon degree-days) for the SNAP2 (mainly residential combustion) emission category, and changed summer/winter ratios for the SNAP-1 (power station) category, a complete revision of the spatial mapping for BVOC emissions (plus an update of the emission factors), addition of soil NO procedures for both global modelling and finer-scale European modelling, addition of seasalt, dust, forest-fires and secondary organic aerosols (SOA) to the standard model. The model has become very flexible, and can now be run with several meteorological drivers, and has the ability to run in nested mode (allowing zooming). A chemical pre-processor has allowed a number of other chemical schemes to be implemented, ranging in complexity from less than a hundred to more than a thousand reactions (the CRI v2 scheme). For sulphur and nitrogen compounds the dry-deposition equations have changed substantially since EMEP-03 (e.g. for the non-stomatal conductance, treatment of humidity, snow, etc.). The aerosol dry deposition scheme is completely new. The MARS equilibrium solver has replaced the earlier EQSAM code, and water associated with PM is now calculated with the same MARS model as used for gas/aerosol partitioning calculations.
Smaller changes include revisions in the equations concerning the stomatal deposition pathway, and in parameter values for land-cover characteristics, the vertical distribution of emissions, collection efficiencies of fine particles, and the VOC speciation. New sources of aircraft emissions and shipping emissions are being used.
In this paper, we have documented the current state of the model, version rv4. The model is continually evolving, but we hope that the rv4 model version provides a good base against which future model changes can be compared. The model code itself is available at www.emep.int, along with the datasets required to run for a full year over Europe.