Acidifying and eutrophying pollutants | EMEP home

The EMEP Lagrangian Acid Deposition Model ( LADM) is a receptor-oriented one layer trajectory model with the spatial resolution of 150x150 km. It includes the chemistry for 10 compounds: NO, NO2, PAN, HNO3, NO3^-, NH4, NO3, NH3 , SO2, SO4^- and [(NH4)2SO4 + NH4 HSO4]/2. Since 1985 the model has been employed at the EMEP/MSC-W to calculate concentrations and depositions of acidifying compounds in Europe, as well as transboundary fluxes and budget matrices.

Content:

The two dimensional mass balance equation for mass concentration q may be seen consisting basically of two parts: one part includes term describing the transport processes, and the other one, Si , accounts for the temporal change of chemical concentration due to all sources and sinks, so that for the component i

A1.1

where u and v represent the horizontal components of wind velocities. In the Lagrangian framework the transport term describes the motion of an air parcel with the wind flow along the prescribed macroscale trajectory. The parcel is assumed to have a height characteristic for the daytime mixing layer and to conserve its volume. Budget equation for the chemistry of an air parcel following the air motions in the boundary layer includes emissions from the underlying grid, chemical depletion and production of the species in the air, and physical removal:

A1.2

here qi is the volume average mass concentration of the component i; Fi represents the concentration production rate, and Ri is the removal term. Production is defined so that

A1.3

here Q is the emission flux density, h is the air parcel height, ki is the chemical production rate of qi from other that i components qj .

The removal term comprises dry and wet removal processes, such that

A1.4

The dry decay coefficient kdi describes surface layer micrometeorology and chemical change:

A1.5

where vdi is the dry deposition velocity and ki is the chemical depletion coefficient. The wet depletion coefficient is expressed as

A1.6

with ( being the scavenging ratio and P the precipitation intensity. It is assumed in the model that during wet periods dry removal occurs concurrently. For a sufficiently small chemical time step (15 minutes step is employed in the model) F and R terms can be considered constants (see Table A1.1).

Calculation of advection, i.e. the trajectory path, follows the methodology described by Pettersen (1956), and that used in the OECD programme on Long Range Transport of Air Pollutants (OECD, 1977, Eliassen 1978). Backward trajectories are calculated at each transport time step on the basis of interpolated wind velocities. From a position, r, at time, t, the initial estimated change in position, Dr0, is defined as

A1.7

The wind vector, v, is interpolated in space between the u and v components in adjacent grid points, and linearly in time between consecutive 6-hourly meteorological input fields. The first guess is corrected as

A1.8

Five iterations are deemed sufficient for convergence to the position r + Dr5. The advection time step length along each trajectory is 2 hours. A total of 49 positions, therefore, are used to define the trajectory path over four days including the start and finish points. Trajectories are defined to arrive four times a day (1200, 1800, 2400 and 0600 UTC) to a set of regularly spaced grid points and to a set of EMEP monitoring stations. In 1997 the calculation domain was extended in the x-direction (see Section A.3.1) so that the entire Turkey, Cyprus and the Mediterranean Sea are covered by the model (Figure A.3.1). Thus, 1260 grid points and 122 measurement sites serve as arrival point in the present LADM version. When the model day is so defined (1200-0600), it corresponds closely to the 24-hour sampling period employed in the EMEP monitoring network, facilitating in this way comparison of the model results with measurement data.

Component	Fi
Nitrogen Monoxide: NO			not wet scavenged
Nitrogen Dioxide: NO2			not wet scavenged
Nitric Acid: HNO3
Peroxyacetyl Nitrate: PAN			not wet scavenged
Particulate Nitrate: NO 3-
Sulphur Dioxide: SO2
Particulate Sulphate: SO 4
Ammonia: NH3
Ammonium Nitrate: NH4 NO3	Non-linear transformations (see A1.4.2 )
Ammonium Sulphate: (NH4)1.5SO4	Non-linear transformations (see A1.4.2 )

A1.3.1 Emission data

Symbol	Definition	Value
QSOx	Emission of sulphur oxides (SOx) per unit area and time	Country-specific seasonal variation
	Fraction of SOx emissions deposited locally in emission grid square	Dependent on emission height and meteorology
	Fraction of SOx emissions emitted directly as particulate sulphate SO4	0.05
QNOx	Emission of nitrogen oxides (NOx) per unit area and time	Country-specific seasonal variation
	Fraction of NOx emissions deposited locally in emission grid square	0.04
	Fraction of NOx emissions emitted directly as NO2	0.05
QNH3	Emission of ammonia NH3 unit area and time	Seasonal sine-variation max=1.3 (June), min=0.70 (January)
omega	Fraction of NH3 emissions deposited locally in emission grid square	w=ww+wd wd=0.19 ww=wd* P/(P+vdNH3/DNH3) wmax = 2 wd

Recent communications with EMEP-CCC imply that this seasonal pattern may not be quite adequate the reality and, therefore, will be a subject of further studies. For emissions from sea areas (international shipping, installations etc.) no seasonal variation is assumed. Marine biogenic sulphur emissions vary with season and latitude with a maximum rate in spring/summer and a minimum in winter. A local deposition correction is applied to emissions as described in Section A1.5.4. Emission notation in Figure A1.1 is defined in Table A1.2.

A1.3.2 Meteorological input

Physical parameter	Level of output	Purpose
From the NWP model:
Horizontal wind velocity (x and y direction)	= 0.925 z 550 m	horizontal transport, eddy diffusivity
Vertical wind velocity	= 0.850 z 1100 m	exchange with free troposphere
Temperature at the ground	z 2 m	aerodynamic resistance
Temperature in the mixed layer	= 0.925 z 550 m	chemical reactions rates
Relative humidity	= 0.925 z 550 m	ammonium nitrate equilibrium, deposition resistance
Cloud cover	free troposphere	photolysis rate of NO2
Precipitation rate	ground level	wet deposition, surface wetness
Surface pressure	= 1. z = 0 m	air density
Turbulent heat flux density	surface layer	aerodynamic resistance
Turbulent stress	surface layer	aerodynamic resistance

Further information on surface properties is required to calculate the dry deposition velocity. Land use data at the resolution of 1/6 x 1/6 degrees was provided by the RIVM (National Institute of Public Health and Environmental Protection, the Netherlands), and subsequently aggregated to the EMEP 50 and 150 km grid at the MSC-W. For dry deposition purposes the classification finally used is divided into 10 different categories: Grass (1), Arable (2), Permanent crops (3), Coniferous forest (4), Deciduous forest (5), Water (6), Urban (7), Other (8), Desert (9) and Ice (10). The two land use classes, "desert" and "ice", were constructed from the NWP-model data (LAM50E). Furthermore, the RIVM class "water" which consists of inland water alone has been revised to include sea areas according to surface area in each grid square. Information concerning snow cover is also employed, with monthly averaged snow cover fields constructed from the NWP-model (LAM50E). The details on the land-use data can be found in Sandnes (1995) and Seland et al. (1995).

Figure A1.1 Overview of the chemical scheme implemented in the EMEP/MSC-W Lagrangian Acid Deposition Model . Q designates the emissions. The calculated components are inside solid boxes. Thin arrows show the chemical pathways; thick solid and dashed arrows are dry and wet deposition respectively. For other notations see Tables A1.2 and A1.4.

A1.4.1 Linear chemical transformations

A1.4.2 Non-linear chemical transformations

and similarly

A1.9

in which [ ]ne indicates non-equilibrium and [ ] e equilibrium concentrations. The production of ammonium nitrate particles, or their evaporation into constituent gases, is determined in the model through estimation of the equilibrium coefficient K:

If non-equilibrium gas concentrations give a product greater than K, ammonium nitrate is created, whilst if the product is less than K the particles evaporate. A quadratic equation defines the equilibrium ammonia concentration, resulting in

  A1.11

with ammonium nitrate and nitric acid concentrations derived from the previous equations. If the ammonium nitrate concentration becomes negative the assumption is of insufficient precursors to reach equilibrium, and remaining constituents are converted to gaseous concentrations. The temperature and humidity dependence of the equilibrium coefficient K is defined according to Stelson and Seinfeld (1982) by comparison between actual humidity and the calculated humidity at the point of deliquescence which marks the point of phase transition. At dry weather conditions, rhd , Stelson and Seinfeld (1982) presented a semiempirical expression for K for solid particulates , in unit of ppm², in a form

  A1.12

When the ambient humidity becomes larger than the relative humidity of deliquescence (rh ³ rhd), Hov et al. (1988) formulated a parameterization of the equilibrium constant for liquid aerosols as

  A1.13

where the relative humidity of deliquescence, rhd , is calculated using an expression derived for the saturated solution by Stelson and Seinfeld (1982) :

  A1.14

The non-linear part of the chemistry asks for a special treatment of country allocation of sulphur and nitrogen that is in the form of ammonium sulphate or ammonium nitrate, whilst for the other components the procedure is straightforward. Is should be stressed that production of ammonium is always calculated from the total concentrations of the precursors and not from individual countries contributions. Country allocation of the particulate ammonium compounds is then performed in two ways: in the first, the allocation is based on the countries' contribution to N from NH3 emissions, while in the second, the origin of S from SO2 emissions and N from NOx emissions determines the allocation.

A1.5.1 Wet removal and sub-grid scale precipitation

.

Symbol	Definition	Value
(i) gases
HNO3	wet scavenging ratio for HNO3	1.4 × 106
HN3	wet scavenging ratio for NH3	1.4 × 10 6
SO2	wet scavenging ratio for SO2	varies seasonaly to reflect variation in H 2O2 ; HN3 = 3 × 105 + 1 × 105 sin[2p (t - t 0)/ta], see note on kt in Table A.3
(ii) particles
NO3	wet scavenging ratio for NO3- and NH4NO 3	1.0 × 106
SO4	scavenging ratio for SO4= and (NH4)1.5 SO 4	1.0 × 106

A1.5.2 Dry deposition

In the LADM, two methods are employed to estimate vd. Reflecting greater scientific understanding of SO2, NH3 , HNO3 and NO2 , the resistance analogy has been introduced for these, with resistances ra , rb and rs parameterised according to sub-grid land use. For PAN and particulates , rb and rs are in effect generalised, and ra calculated accordingly. Application of the resistance method in the model is described in detail in Jakobsen et al., (1996).

Aerodynamic Resistances, ra .

Aerodynamic resistance, ra, is described in the standard way ( Garland, 1978) by

 A1.15

Here z0 is the roughness length; L is the Monin-Obukhov length defined as standard (e.g. Stull, 1988), using values for turbulent stress, t, heat flux density, surface air temperature and pressure, r, from the Numerical Weather Prediction model; k is von Karman's constant using a value of 0.37, and u* = (t/r)^0.5 is the grid averaged friction velocity, and d is the displacement height taken as 70% of vegetation height. The roughness heights are assumed to be 10% of vegetation height (Table A1.7). For PAN and particulates the sub-grid surface resistance is not parametrised explicitly and ra given by Equation A1.15 is used. For the others (SO2, NH3, HNO3 and NO2), roughness length varies with sub-grid land use, and hence typical friction velocity must be redefined. From the NWP model only grid square averaged values of u*, T*, L and z0 are available. These parameters are then converted to land use type specific values u*lu, T*lu and Llu within each grid square applying the land use type specific values of z0,lu from RIVM (Table A1.7). Finally, the land use specific aerodynamic resistance ra,lu is derived :

 A1.16

Here Y is the stability function accounting for diabatic corrections to the logarithmic profiles of momentum and heat. The algorithm converting the grid square average values to land use type specific values within each grid square is further described and discussed in Jakobsen et al., (1996).

Over water surfaces the roughness length is given by Charnock relation :

 A1.17

The value of Charnocks constant is taken from the NWP model as described by Nordeng (1986) and Nordeng (1991).

SO2, HNO3, NH3 and NO2 quasi-laminar boundary layer resistance, rb .

Over land surfaces, the resistance of the quasi-laminar layer to transfer, rb, is determined by scaling of the friction velocity according to the diffusivity characteristics of the gas. Following Hicks et al. (1987), the resistance above land is given as:

 A1.18

in which Sc is the Schmidt number (defined as n/Di with n being the kinematic viscosity of air (0.15 cm² s-1) and Di the molecular diffusivity of gas i), and Pr=0.72 is the Prandtl number.

Over water surfaces (i.e. sea) the quasi-laminar resistance r b is approximated by the procedure presented by Hicks and Liss ( 1976) :

 A1.19

The problem that rb may become as negative that very large and also negative deposition velocities appear has to be prevented. Based on physical considerations, the limitation 0 ms^-1 < vd < 0.1 ms^-1 is imposed to avoid occurrence of large and negative deposition velocities..

SO2, HNO3, NH3 and NO2 Surface Resistances , rs .

The resistance to uptake or destruction at the surface, rs , depends on a combination of the biological and physico-chemical properties of the absorbing surface (such as pH, stomata, surface characteristics, etc. ) and the properties of the component. The surface resistances used here are mainly derived from Erisman et. al., (1994).

For SO2, NH3 and NO 2 on vegetation not covered with snow the surface resistance can be expressed as

For non-vegetative surfaces and any surfaces covered by ice and snow r s = rsoil and prescribed in Table A1.11.

For HNO3 resistance is considered essentially aerodynamic, and rs is set to 10s m ^-1 ( 50s m^-1 when both -5 ^oC and snow covered).

For non-vegetative surfaces, rs is more simple and the values are given in Table A1.11.

Remaining Chemical Species (PAN and particulates)

For PAN and particulates for which detailed resistance data is less available a dry deposition velocity comparing with reported measured values is determined for 1 m. First, the maximum deposition velocity, vdmax, is selected (Table A1.12). Then the vdmax is adjusted in a regular way to reflect seasonal surface affinity change such as observed reductions in deposition rates over snow and expected higher vegetation uptake in summer (e.g. Whelpdale and Shaw, 1974). This seasonal/latitude adjustment is applied so that the actual dry deposition velocity at 1 m, vd1 < vdmax on land through

  A1.21

Here B represents latitude, t is the time of year, and t is the time of day. The latitudinal adjustment is produced by the function f(B) = r/R where r is the distance from current position to the North Pole and R is the distance between equator and North Pole. The multipliers a1 and a2 are the trigonometric functions

in which t0 is February 1st, and ta represents one year. The daily variation of vegetation uptake is applied simply with the multiplier D set to 1 between 0400 and 2200 local time, and to 0.25 for the remaining 6 hours. No variations are applied to sea deposition, although the current formulation better reflects summer conditions when relative contributions from dry deposition may be at their most important (Barrett, 1994). Particulates dry deposition is assumed to be much slower, given typically small particle sizes. A constant of 0.1 cm s^-1 is applied to sulphate and nitrate aerosols, approximating a size of 0.1-1 mm (e.g. Sehmel, 1980).

Dry depletion is applied above the surface layer, and so the appropriate 1 m value is transformed to a 50 m (typical surface layer depth) value by application of similarity theory, thereby reflecting the effect of aerodynamic resistance between 1 m and 50 m (Eliassen and Saltbones, 1983):

 A1.23

No profile adjustment is made to the vd1 set for aerosols, i.e. vd50 = vd1 .

The determined deposition velocity assigned to a grid square is an area average, reflecting the percentage occurence of various sub-grid landuse types within the square. Hence:

 A1.24

where vd,n(wet) and vd,n(dry) are dry deposition velocities for wet and dry parts respectively of land use n; an is the fraction of land use n in the grid square; and ( is the fraction of wet part.

A1.5.3 Combined effect of wet and dry depletion

 A1.25

A1.5.4 Local deposition of emissions.

The assumption of instantaneous mixing of emissions throughout the air parcel is not realistic. If no correction is made, the model will very often underestimate near-ground concentrations before complete mixing is achieved, particularly close to emission sources, and with this will underestimate the dry deposition flux. To compensate that a fraction of emissions is assumed to dry deposit directly within the emitting grid square in addition to what is calculated in a normal way, with the remainder available for mixing, transport and further depletion. The height of emission source is critical to the local deposition fraction. An approximation is necessitated, by which the proportion of NH3 emissions additionally locally deposited is taken as 0.19, and of NO2 as 0.04. As gradients from ground level sources of ammonia emissions may be very strong, and NH3 is highly soluble, a correction on wet local deposition is also made to NH3 on the basis of the scavenging ratio and precipitation (Iversen et al., 1990), to a maximum of a factor of 0.19. As a first guide to current model estimates, Janssen and Asman (1988) have suggested that similar values under neutral conditions may be appropriate for NH3 and NO2 emissions from heights of 0 m and 100 m respectively. A more detailed approximation is made for SO2 derived from the work of Tuovinen and Krüger (1994).

The local deposition factor is calculated for each grid square on the basis of vertical emission distribution and meteorology. With a 'low/high' division of emission sources defined as 100m, the local deposition factor is determined separately for each fraction, so that the overall factor, atotal, may be described by:

 A1.26

The derivation of the separate a factors is taken from look-up tables quoted in Lehman (1991). These were calculated assuming a 1m deposition velocity of 8 mm s^-1, and so must be scaled according to actual conditions, thus:

 A1.27

in which vdREF represents a scaling factor. Starting from the deposition velocities already calculated for SO2 and equation (17) in which vd(1m) is set to 0.008m.s^-1:

 A1.28

Higher 1m deposition velocities than the 0.008m s^-1 assumed by Lehman (1991) are now possible, so that the a LEHMAN values are extrapolated when necessary, to an artificial six-hour upper limit of 0.50. Finally, emission data available to the MSC-W is as annual totals, which is then scaled with a monthly cycle. The a values are, therefore, calculated as monthly mean high and low values from the six-hour meteorology. The local deposition defined from (23) is allocated in each square according to the emission country, i.e. each of the countries in a square receives the whole additional deposition originating from emissions from that country within the grid.. This procedure affects only the pollutant budgets, not the geographical fields.

The EMEP model is a single-layer model representing the mixed layer over Europe and the eastern North Atlantic. However, there is no physical isolation of this layer from the remaining atmosphere, therefore, mass exchange both with the free troposphere and areas beyond the model domain should be accounted for. Regarding boundary layer meteorological parameters, the EMEP model is a sub-domain of the Numerical Weather Prediction Model of DNMI, so that meteorological estimates, especially over sea areas where observations are scarce, are fully developed before encountering the EMEP domain. For chemical characteristics, three types of data must be pre-determined. These are: a) free tropospheric conditions, b) background concentrations of non-calculated components, and c) initial concentrations for calculated components. Since the origins of these concentrations cannot be traced by the model, they are attributed to concentrations and depositions of indeterminate origins.

A1.6.1 Exchange with the free troposphere

Even though most transport in the model is expected within the mixed layer, intermixing with the tropospheric concentrations above the modelled layer may be of a significant importance, particularly under transport to long distances. Therefore, the interaction with the free troposphere is included necessarily to the LADM as described by Eliassen and Saltbones (1983). As the height of the air column changes it is compared with the analysed maximal mixing layer height: if the box is higher than the local mixing height, the upper part of the box is assumed to be injected into the free troposphere and is not followed further in the model; if the opposite is true, tropospheric background concentrations are mixed down into box. The free troposphere is approximated, thus, by an infinite reservoir with constant low concentrations.

Mixing heights are calculated from the observed temperature vertical profiles. Estimated thus mixing heights across Europe at 1200 UTC are interpolated across the domain to define the initial mixed layer height at t1. The height of the advecting air parcel is assumed to change along its trajectory according to vertical motion due to horizontal divergence/convergence so that its volume remains conserved. At t2 = t1 + 24 h the height is

 A1.29

where wr is the vertical velocity on top of the mixed layer, interpolated in time from the 6-hourly input fields, and reduced linearly for heights below 1000 m. At t2 the height of the air parcel is redefined to equal the local mixing height Hm(r(t2), t2). The difference (Hm- h) represents an exchange volume, so that the revised parcel concentrations q' become

 A1.30

in which qa is the concentration in the free troposphere.

Contributions arising from free troposphere exchange contribute in the model to the concentrations of indeterminate origin. The total indeterminate contribution to deposition ranges between 5 and >20% (Sandnes, 1993), of which about 80% is estimated to come from the free troposphere.

A1.6.2 Initial and boundary concentrations for calculated components

The chemical characteristics of the boundary layer air mass must be described at the start of the trajectory. For start positions within the calculation domain (36x35 points), the initial values can be obtained either from the previously calculated fields through assimilation of these values, or by prescribing them in case of initialising the start of a run period. Outside the calculation domain, boundary layer and free troposphere concentrations are to be given at any time. The initial and boundary concentrations of sulphur compounds (SO2 and SO4) are given based on calculations from the 3-D hemispheric model (Iversen and Tarrason, 1990) with adjustment to reflect measurements from a variety of sources, although such measurements are generally sparse and discrete. For the sum of the three nitrate forms, the total number density within the boundary layer is assumed to equal that for total sulphate. reasoned by the fact that both are secondary species originating from emission sources of similar patterns. However, nitrate is deposited faster than sulphate, such that with ageing nitrate concentrations should be lower than those of sulphate. Therefore, nitrate concentrations in the free troposphere is assumed to equal half the total sulphate concentrations. The speciation of nitrate and sulphate is determined by immediate turnover and equilibrium as described in A1.3.2. Furthermore, the review of measurement data by Fehsenfeld et al. (1988) for boundary layer NOx and PAN has been used to assign the initial and boundary concentrations. Finally, for ammonia the concentration pattern earlier estimated by the LADM is adopted as a guideline. It is used together with the assumption that ammonia concentrations must closely reflect emission distribution due to its rapid depletion. Above the mixing height concentrations are assumed to be 10% of the boundary layer values. The standard chemical scheme provides the speciation. As for tropospheric concentrations, Sandnes (1993) has previously illustrated some so-called background values used in this procedure. In general the relationship between concentrations and temporal emission pattern is much stronger within the boundary layer for both nitrogen and sulphur compounds. These background values were estimated on a monthly basis for generalised regions of similar pollutant loading. For sulphur compounds Central Europe, Fennoscandia, North Africa, remote European Commonwealth of Independent States, and the North Atlantic Ocean are so defined, whilst with respect to NOx where the estimates are less certain only three regions are defined: Central Europe, remote continental Europe & North Africa, and the North Atlantic. Through the year the NOx concentrations mirror the assumed general temporal variations in emissions, whilst for SO2 and total sulphate very little seasonal variation is evident.

A1.6.3 Background concentrations of non-calculated components

Three components which are not calculated in the model, but required to estimate reaction speeds are taken from prescribed fields, these being ozone O3, and the hydroxyl OH and peroxyacetyl CH3COO2 radicals. The values have originally been obtained from the 2-D global model of Isaksen and Hov (1987). Comparison has then been made with reviewed measurement data which indicates some model overprediction, in particular winter and spring ozone levels at high altitude mountain sites representative of the background atmosphere (Logan, 1985). A latitude dependent correction has, therefore, been applied to the seasonally varying ozone concentrations in order to ensure better coincidence with the observations. Additionally, the minimum ozone mixing ratio of 25 ppb was imposed. Furthermore, the two radicals have been adjusted with a diurnal cycle parallel with sunlight. Some examples of values are given in Table A7 below.

References