(2.1)The state of the troposphere is determined by surface forcings. Weather systems are driven by heating and moisture release from the surface and act throughout the depth of the troposphere. However, only a relatively thin layer of this lowest part of the atmosphere is regularly subject to turbulent mixing, the subject addressed here using micrometeorological methods. This layer is known as the boundary layer. The defining features of the boundary layer are that it responds to surface forcings on a time scale of around an hour, and it is turbulent throughout its depth due to interactions with the surface (in fact these two are linked because the turbulence allows rapid equilibration with the surface). The following sections introduce some fundamental concepts for study of transport in the boundary layer.
2.1.1 Turbulence Structure
Boundary layer turbulence is generated on a variety of scales. At night the scale of turbulence generated by wind shear in the nocturnal boundary layer may be as small as 100 m. During a summer’s day thermally generated turbulent motions can span the whole depth of the mixed layer (up to 4000 m in mid-latitudes). In between these scales are wake–generated turbulence and smaller scale thermal motions. At smaller sizes the combined effects of wind shear and shearing of large scale turbulence produce smaller turbulent motions, which in turn feed even smaller eddies. At scales where the viscosity of gases becomes important, turbulent motion dissipates as heat.
The process whereby energy is put into the largest scales of turbulence, passed through smaller and smaller eddies and finally dissipated into heat is known as the energy cascade. In this chapter, reference will be made to the primary (or energy containing) motions and also to the behaviour of the smaller eddies which have no source or sink of energy except eddies larger or smaller than themselves. This range of turbulence scales is known as the inertial sub–range and can be probed using Fourier analysis.
2.1.2 Boundary Layer wind speed profile
Wind speed varies throughout the boundary layer depth. Near the surface the wind speed is zero, above the boundary layer it is approximately geostrophic, and the transition is approximately logarithmic with height. The change in speed is due to vertical mixing of momentum from aloft down towards the surface. The mixing is caused by turbulent motions which have characteristic velocities of the order of (for the non-advective part) u’, v’ and w’ (this notation is explained in the section on eddy covariance) in the x, y and z directions respectively. These can be approximated by:
(2.1)
Where l = mixing length (m), z =
height above the surface (m), U = mean wind speed (m s-1),
but note that later
is also used for consistency with eddy covariance nomenclature and u*
= friction velocity, the tangential velocity of the eddy motions (m
s-1). The mixing length is the characteristic size of
eddies in the surface layer, and the friction velocity can be
interpreted as the mean wind speed change over one mixing length. The
change in speed with height gives rise to a Reynolds stress as
parcels of air with finite vertical extent are deformed. The
magnitude of this (one-dimensional) stress can be calculated from:
(2.2)
Where = Reynolds stress (kg m-2 s-1) and = density of air (kg m-3). The mixing length, l is a strong function of height and may be approximated by:
(2.3)
Here k = the von Karman constant, d = zero plane displacement (m) and m = the empirical stability correction (discussed later). Figure 1 shows schematically the relationship between wind speed profiles, mixing length and stability (discussed later). The zero plane displacement describes the increase in apparent ground level over a rough surface. It is the theoretical sink level for momentum, where the momentum flux reaches zero. It has a non–zero value over all but perfectly smooth surfaces, and is significantly higher for forest than for grassland. Thus (z - d) represents height above ‘aerodynamic ground level’. The von Karman constant is used often in micrometeorology and has an empirically determined value of 0.41 (dimensionless). It can be interpreted here as the mean eddy size per metre above ground level.
It is worth noting here that since the mixing length increases with height (figure 1), so does the time–scale on which “primary” turbulence operates. This becomes significant with regard to eddy covariance. According to Taylor’s (‘frozen turbulence’) hypothesis and from (2.3), if Te is a turbulence time–scale, (=l) a length scale and U the wind speed:
(2.4)
Returning to the wind speed gradient, combining (2.1) and (2.3) gives:
(2.5)
Which can be solved for dU and integrated over a height range treating m as constant with height to yield the logarithmic wind profile:
(2.6)
Where the constant of integration that has appeared (z0) is the height above d where the windspeed theoretically reaches zero. It is called the aerodynamic roughness length and has units of metres. It is a difficult quantity to measure, as it is sensitive to the value of d used as well as m, not the most straightforwardly applicable atmospheric variable for reasons which will be discussed. At this point it is worth noting that the logarithmic wind profile is purely a result of the Reynolds stress or momentum flux. The wind speed is zero at the surface (as with all boundary layer flows) and momentum is mixed down through the resulting gradient from aloft. The stability correction is a result of the fact that the rate of vertical momentum mixing (see later section on eddy diffusivity for momentum) is a function of atmospheric stability, as shown in figure 1. The cause and effects of changes in stability are examined in the following section.
2.2 Atmospheric stability
2.2.1 Introduction
Due to the adiabatic expansion of a parcel of air as it rises, air raised to higher altitudes has a lower temperature than air of the same composition and initial temperature at lower altitudes. This results in a fall of air temperature with height called, for dry air, the dry adiabatic lapse rate. This lapse rate is 9.76 K km-1 (determined using the ideal gas equation and an assumption of hydrostatic balance). Due to surface forcings, the ambient lapse rate is rarely adiabatic. Where the ambient temperature lapses faster with height than the adiabatic lapse rate, an air parcel moving upwards adiabatically becomes warmer and hence less dense than the air around it. This causes it to continue rising more rapidly until reaching an inversion (change to a lower lapse rate) in the ambient temperature gradient. An atmosphere where the lapse rate is super–adiabatic is said to be unstable. Conversely if the ambient lapse rate is sub–adiabatic, a displaced air parcel tends to return to its original height and vertical mixing is suppressed. This is the stable case. Where the ambient lapse rate is exactly adiabatic, the atmosphere is said to be neutral, and the vertical mixing caused by mechanical turbulence is neither thermally suppressed nor enhanced. Figure 1 demonstrates the effect of stability on the wind speed profile due to changes in mixing length.
2.2.2 Virtual Potential Temperature
Due to the lapse of temperature with height, it is convenient to define a variable such as a height independent temperature. This is again done using the ideal gas relation such that:
(2.7)
(2.8)
Where Ta = air temperature (K), p and p0 = pressure at the measurement height and reference height respectively (mb), = potential temperature (K), cp = specific heat capacity of dry air at constant pressure (1004.67 J kg-1 K-1) and cv = specific heat capacity of dry air at constant volume (717.02 J kg-1 K-1). Conventionally p0 is taken to be 1013.25 mb (standard atmospheric pressure) or the pressure at the surface (when this is measured). Potential temperature can then be interpreted as the temperature an air parcel would have when moved adiabatically to the standard atmospheric pressure level or to the surface respectively. This allows easy assessment of the thermal energy of an air parcel.
Potential temperature is also a useful measure of the potential for buoyancy of air. If dry air has a higher potential temperature than its surroundings it will be more buoyant. This can be generalised to moist air by defining the virtual potential temperature. Virtual temperature is the temperature air must have to be at the same density as moist air at the same pressure. Combining these two concepts, virtual potential temperature is defined (for saturated air) as:
(2.9)
Where rsat = is the saturation mass mixing ratio (kg kg-1, g g-1) water in the air and rL is the actual water vapour mixing ratio. Using this definition air at any level and with any water content can be compared for buoyancy with other air.
2.2.3 Obukhov Length and the Stability Parameter
The stability at a certain height could be gauged by the virtual potential temperature profile at that height, but this is not a convenient quantity to measure and would not give information on the stability of the whole depth of the boundary layer. Normally as a measure of stability we use the Obukhov length, which is “…the characteristic [height] scale of the sub–layer of dynamic turbulence” (Obukhov, 1946). This interpretation also means the height up to which mechanically generated turbulence dominates over thermally induced turbulence. The Obukhov length (in m) is defined as:
(2.10)
Where g = acceleration due to gravity (9.81 m s-2) and H = sensible heat flux (W m-2). The Obukhov length is the most important scale for stability in the surface layer, and is used widely in this work. It is worth noting that the virtual potential temperature used here is not derived from equation (2.9). The results presented were all obtained using sonic anemometry, and since the anemometer measures air density and converts this to a temperature, this direct determination of virtual temperature was used.
The Obukhov length takes values from around -2000 to 2000 m. Equation (2.10) makes it clear that in stable conditions (negative heat flux) L is positive and in unstable conditions it becomes negative. Strong winds producing large friction velocities tend to increase its value. Large moduli of Obukhov length indicate near neutral conditions. In the neutral limit the value of L is infinity (division by zero heat flux in equation 2.10). This gives rise to strange behaviour in time series of L. Instead of using L to describe stability this work adopts the convention of using the stability parameter (n.b. NOT the same as the empirical stability correction defined in (2.12) and (2.13)). The stability parameter is defined as:
(2.11)
The stability parameter is dimensionless. It has a value of zero in the neutral limit, is positive in unstable conditions and becomes negative in stable conditions.
2.2.4 Monin - Obukhov Similarity Theory
Monin - Obukhov similarity theory is the basis of much of the micrometeorology used here. It applies to the surface layer (approximately the lowest 10% of the boundary layer) and is often referred to as surface layer similarity. In situations where the mechanisms determining the behaviour of an atmospheric parameter are incompletely understood, similarity can be used to derive empirical relationships, typically between height and the parameter in question. The approach is to non–dimensionalise groups of variables using scales such as the Obukhov length, and to look for common (“similar”) behaviour in the non–dimensionalised variables. Other parameters often used to normalise variables include wind speed, measurement height, friction velocity and roughness length. A similarity method is used in the development of the empirical model for urban fine aerosol concentration, which is presented in chapter 5.
2.2.5 Townsend’s Hypothesis and Similarity
Townsend’s hypothesis states that boundary layer turbulence can be divided into an active, momentum transporting part and an inactive part. It further states that these two forms of turbulence do not interact. Active turbulence is created by wind shear and has properties scaling on local (surface layer) variables. The “inactive” turbulence is the result of the energy cascade, – removal of energy from the active structures and redistribution down to dissipative scales. The processes leading to the energy cascade are remote from the surface and are thought to scale with so called “outer layer” parameters (McNaughton and Brunet, 2001) such as the boundary layer depth.
The distinction between these two “types” of turbulence gives rise to a potential problem with Monin – Obukhov similarity theory. Active turbulence (the energy containing part) should scale well with the surface layer parameters listed in the previous section, but secondary or inactive turbulence might not be well represented by surface layer scaling. In fact, McNaughton and Brunet (2001) state that Monin – Obukhov similarity appears only to be fully successful under turbulent stable conditions (e.g. the nocturnal boundary layer) where eddy covariance measurements of the type presented here tend to be least accurate anyway, due to the low magnitudes of fluxes observed. However, Monin – Obukhov scaling is used in this work in the absence of a more reliable framework.
2.2.6 Stability Correction
As noted previously (figure 1), stability affects the rate of exchange of momentum, altering the wind profile. In fact it affects vertical transport of all entrained properties. The empirical stability correction m was introduced as a way of correcting predictions of wind profiles for non-neutrality. There are also corrections for temperature gradients and trace substance profiles. These are H and respectively. They are equal to one another in stable conditions, while for neutral conditions the momentum correction is different. The values of have been the subject of much investigation. The values used here are those originally obtained by Dyer and Hicks (1970) and Businger (1966) for unstable conditions and from Webb (1970) for stable conditions.
L
< 0 (2.12)
L
> 0 (2.13)
In fact even with the use of these corrections the profiles of wind, temperature and trace concentration are not perfectly logarithmic. Strictly the stability corrections are themselves a function of height, and m should be integrated with height in the derivation of (2.5). However, this complicates the analysis significantly (see e.g. Panofsky, 1963), and has not been applied here because the use of the wind profile is limited to calculations of z0 and d. Since most of the data used for these calculations were taken in near stable conditions the error was allowed in order to simplify the analysis. The theory outlined here is known as the approximate log–linear approach, and works well for point measurements in conditions of “non–extreme” non–neutrality.
Having gained an impression of the structure and basic dynamics of the boundary layer the following sections outline some of the micrometeorological techniques applied in the studies presented here. The following sections outline the gradient and eddy covariance methods. Later, limitations of the techniques are discussed and further analysis methods considered.
2.3 Gradient Flux Measurements
2.3.1 Introduction
Measurements of scalar fluxes can be made using information on the vertical profile of the scalar and an estimate of the eddy diffusivity for that quantity. Eddy diffusivity is merely a conceptual tool, but allows estimation of fluxes of momentum, heat and entrained quantities from wind speed and concentration profile measurements. It is based on the concept of turbulent diffusion of a property along its mean gradient. Turbulent diffusion is used in analogy with molecular diffusion, except that where turbulence exists in the atmosphere it is orders of magnitude more efficient at transporting properties than molecular diffusion. This method has advantages over inferential methods such as chamber or cuvette measurements since it represents a flux due to fully developed turbulence over a large spatial scale (depending on height, approximately of the order of hectares). The following sections give a brief outline of the technique, but do not go into great depth because although gradient fluxes are presented in this work, they are used only to support and interpret results from other techniques, and were not an integral part of the work.
2.3.2 Gradient Method Theory
The basic equations used to calculate fluxes using the gradient method are presented below. They are all based on the eddy diffusivity for the quantity of interest multiplied by the vertical profile of that quantity.
(2.14)
(2.15)
(2.16)
(2.17)
Where H = sensible heat flux (W m-2), E = latent heat flux (W m-2), f = entrained species flux (e.g. aerosol is in cm-2 s-1), = density of dry air (1.225 kg m-3 at sea level), cp = specific heat capacity of dry air at constant pressure (1004.67 J kg-1 oC-1), and Km, KH, KE and K are the eddy diffusivities for momentum, heat, water vapour and entrained properties respectively. The gradients can be estimated from measurements at several heights, and all other factors in equations (2.14 – 2.17) are constants. This just leaves the values of eddy diffusivity to be calculated in order to determine fluxes.
There is evidence (Thom, 1975) that in all stability conditions, KH = KE = K, and that in neutral and stable conditions these are equal to Km. In unstable conditions, they are not equal to Km due to the difference in values of the stability correction. Using (2.14) with (2.1 – 2.3) it can be shown that:
(2.18)
This relationship can be used to estimate all the eddy diffusivities in stable and neutral conditions, but for unstable conditions an analogy must be made using the definitions of m and H so that:
(2.19)
Note that the friction velocity in (2.19) can still be found using (2.18) in unstable conditions. This formulation of the gradient technique has only one remaining obstacle. The values of m and H as defined in (2.12 – 2.13) contain the Obukhov length, and according to (2.10) this requires the values of H and u* as input. A relationship exists between easily measurable quantities and the stability correction, and this can be obtained by replacing the term containing L again from Thom (1975):
(2.20)
Where Ri is the gradient Richardson number. It is known as the gradient Richardson number because the definition used here is in terms of gradients of temperature and wind speed:
(2.21)
Here T = air temperature (oC) which is normally measured at similar heights to the scalar and wind speed, and g = acceleration due to gravity (9.81 m s-2). It is clear from the level of complexity of the defining equations that the gradient method is not trivially easy to employ. Some researchers use the technique for its ability to derive fluxes from slow scalar sensors, but also employ sonic anemometry to directly measure H and u*, making the last two relationships (2.20 – 2.21) redundant and greatly simplifying the analysis.
The definitions of the gradient theory presented here are, in fact, the simplest possible implementation. A more exhaustive survey including methods for linearising the measured profiles with height and a more accurate treatment of the stability correction (which is in reality height dependent) is presented in Sutton (1990).
2.4 Eddy Covariance
2.4.1 Introduction
The defining feature of the boundary layer is that it is turbulent throughout its depth. This means that entrained entities such as heat, moisture, pollutants and momentum are transported very efficiently in the vertical. This turbulent transport dominates other transport mechanisms such as molecular diffusion for all but the largest aerosol, which can sediment out rather rapidly. It is therefore useful to be able to quantify turbulent transport. There are two widely used micrometeorological techniques for flux measurement, although other inferential techniques exist (e.g. chamber studies and box or mass balance methods). A third micrometeorological approach, relaxed eddy accumulation exists. However it is not used in this work and is treated adequately in Gallagher et. al. (2000). The two common direct methods are the eddy covariance technique and the aerodynamic gradient method.
The major advantage of the gradient method is that since scalar concentration does not have to be measured rapidly enough to distinguish the effects of individual eddies, slow response sensors can be used. In this work, fast enough sensors were either available or were developed to use the eddy covariance technique for aerosol fluxes. This approach was followed for two reasons. First, eddy covariance measurements are made at a single height, meaning that only turbulent transport contributes to the measured flux. In the gradient technique, any divergence in the flux (e.g. chemical reactions of ammonia (chapter 4) or aerosol coagulation (chapter 5)) over the height interval of the gradient measurement would result in a gradient not representative of transport of the species and hence a spurious estimate of the flux. The second advantage of eddy covariance is that since the measurements are made at high frequency and resolve the effects of individual eddies, spectral, co–spectral and probability density analyses (see later sections) can be employed to determine the mechanisms responsible for the flux.
For eddy covariance the three–dimensional wind field is measured rapidly, and the value of the scalar of interest (e.g. temperature, water vapour mixing ratio, aerosol concentration) is measured at the same frequency and location. The covariance of the vertical wind speed and the scalar value is equal to the flux at the measurement point. Covariances are calculated over periods of fifteen to thirty minutes depending on stationarity issues (discussed later). The following sections define the eddy covariance flux and discuss some requirements and limitations.
2.4.2 Eddy Covariance Definition
The covariance of two quantities is defined as:
` (2.16)
Where rxy = covariance of x and y and the overbar denotes an average. N = number of data in the time series, in this work N was determined by the logging frequency and the fifteen minute averaging period always used. The time series are subject to Reynolds decomposition before the covariance calculation such that:
and 
(2.17)
Where
= average of time series x over the period of the calculation
and x’ is known as the perturbation part of x.
Similar definitions for y apply. Examples follow to illustrate
the use of eddy correlation to calculate sensible heat flux and
aerosol flux.
2.4.3 Sensible Heat Flux
The sonic anemometer (see chapter 3 on the CPC flux system) calculates the virtual potential temperature at the same time as the wind speed using a measurement of the speed of sound. In this way the measurements are coincident and co–located. Replacing x with the vertical wind speed and y with the temperature data gives:
(2.18)
Where rw = kinematic sensible heat flux, more commonly known as H (W m-2), = density of dry air (kg m-3) and cp = specific heat capacity of dry air at constant pressure. Note that the constants are required to convert the kinematic heat flux into the required units. The units of covariance without modification are the rather less intuitive oC m s-1. The overbar notation is used from now on as it is the most widely used in the literature. However, for computational efficiency the following relation was used wherever the overbar covariance notation appears.
(2.19)
The eddy covariance technique works by averaging the effect of many eddies over time. In this example, the movement of heat is detected as follows. If the temperature gradient is such that turbulence tends to move heat upwards, when air is moving upward it tends to be warmer than average. Descending air tends to be cooler than average, and the product inside the summation in equation (2.16) is most often positive, giving a positive value for H. Where heat is being transported downward the opposite is true.
2.4.4 Latent Heat Flux
Measurements of latent heat flux are conducted in a similar way to sensible heat flux. A separate instrument is co–located with the sonic anemometer (within around 10 cm) to ensure that both instruments are measuring properties of the same turbulent structure. This instrument, a Krypton hygrometer (KH2O, Campbell Scientific Inc., Logan, Utah, USA) which uses extinction of UV light to measure fluctuations in vapour concentration (Campbell and Tanner, 1985), outputs data in mV at the same rate as the anemometer. The calculation then becomes:
(2.20)
Where E
= latent heat flux (W m-2), Lv is the
latent heat of vaporisation for water and q is the water
vapour loading (kg m-3). The Krypton Hygrometer cannot
measure the absolute value of q because the conversion for the
q signal (Vq) from volts to (g m-3)
is subject to significant drift due to fouling of windows in the
optical path. Since the shape of the calibration (
)
is very stable, the measurements can be used to directly determine
Vq’. Equation (2.20) then becomes:
(2.21)
With Vq in mV, x =
optical path length of the instrument (cm),
= average signal (mV) and kw = extinction
coefficient (mV m3 g-1 cm-1). In
fact, the extinction coefficient kw is determined
three times during the manufacturer’s calibration. The first
value provided is an average for the full range of water vapour
density (dry to saturated), while the other two values are used for
the “wet” (8–20 g m-3) and “dry”
(1.5–10 g m-3) ranges. The separation increases the
accuracy of the measurement by allowing for a slight
non–linearity in the response of the instrument. In this work,
a separate humidity sensor (Humicap HMP 50Y, Vaisala, Helsinki,
Finland) was used to determine which of the two latter calibration
ranges to use.
2.4.5 Aerosol Number Flux
Aerosol number flux calculations, both total fine mode and size segregated accumulation mode, are outlined in more detail in chapters three and four respectively. The basic approach (for total number flux) is to calculate the covariance of vertical wind speed and aerosol number concentration, in a manner analogous to the heat fluxes.
(2.22)
Where F = aerosol number flux (cm-2 s-1), w = vertical wind speed (m s-1) and = aerosol number concentration (cm-3). The application and interpretation of eddy covariance aerosol number flux measurements are discussed later in this chapter and in the results sections. The operational aspects of all types of eddy covariance measurements conducted as part of this work are covered in the later chapter on the aerosol flux measurement system.
2.4.6 Vertical Sensor alignment
Equation (2.17) shows the Reynolds decomposition used in eddy covariance calculation. For the calculation of the perturbation part of the w series (w’) it is essential that the z axis of the anemometer be normal to the mean wind. If this is not the case, gusts in the x and y directions (u’ and v’ fluctuations) contaminate the w’ series. The calculated flux is then representative of turbulent transport not in the vertical, but in some ill–defined direction that includes a component of the longitudinal (advective) and transverse flows.
In the field it is near impossible to align an
instrument perfectly vertically, and in any case the mean flow is not
horizontal except over perfectly flat terrain. The solution is
to rotate the co–ordinate system of the wind speed measurements
so that by definition
.
In this newly defined co–ordinate system the measurement
represents the turbulent flux perpendicular to the mean flow. This is
achieved using the following, where the subscript “o”
denotes the rotated quantities that satisfy the condition above.
(2.23a)
(2.23b)
(2.23c)
The angles and are defined by the following, where the over–bar, as usual represents an average:
(2.24a)
(2.24b)
In this work, calculations (2.23 a–c) and (2.24 a-b) were used and the Reynolds decomposition performed on the resulting data series (uo, vo, wo) for simplicity. Another method exists in which fluxes can be calculated and the rotations performed on these single values rather than repeatedly over the course of a whole time series. This method has been used widely in the past due to its computational efficiency, but was not used here for the sake of simplicity. For completeness the algorithm is presented below, with the same definitions for and as in (2.24).
(2.25)
2.4.7 Sensor Requirements
For eddy covariance to work there are certain requirements for the sensors to be used. These are primarily based on sensor speed. In order to measure a turbulent flux, the sensor must respond to changes at a rate fast enough to resolve the flux. This means that if a property is being transported by eddies with a time–scale of (e.g.) ¼ s, the sensor must have a response time better than ¼ s. This is not a problem for sonic anemometers at the heights used in this work, but can be a problem within canopies (recall that the mixing length and hence mean eddy size increase with height – equation 2.4). A more serious problem is encountered with chemical and aerosol sensors. These can have much longer response times (see chapter 3) and for some species sensors do not yet exist which are fast enough to perform eddy covariance (e.g. ammonia). More detailed response requirements are outlined in the section on spectral data analysis.
There are several methods of correcting calculated fluxes according to the frequency response of the instruments used. Eugster and Senn (1995) develop a model in which spectral and co–spectral analyses are used in conjunction with an analogy with electrical inductance to correct flux measurements for damping in tubes and analysers. In another approach, Massman (2000) uses co–spectral transfer functions to describe and correct for the effects of damping, sensor response and even the averaging period used. A discussion and subsequent revision of this paper were published by Rannick (2001) and Massman (2001). The necessity for the use of such techniques is assessed on an individual basis for each data set presented in this thesis, - see individual results sections for details (chapters 4 and 5).
Another requirement for eddy covariance sensors is that they must have sufficient resolution to measure the expected changes in concentration. The expected variation may be estimated with knowledge of the gradient throughout the boundary layer, or from standard deviations reported from other studies. For example, it would be pointless to attempt to measure a gas flux resulting from a 50 ppb gradient over the mixing length using a sensor with 1 ppm resolution. The turbulent fluctuations in the concentration of an entity may be extremely small, but over a period still amount to an environmentally significant flux. This leads on to the matter of counting statistics. This is not so much of a problem in measuring gas fluxes, but the discrete nature of aerosol means that sufficient individual aerosol particles must be observed over the period of a measurement so that the ‘one over root N’ error (assuming Poisson counting statistics) does not become significant. Fairall (1984) estimates that the minimum number of observed aerosol required for a statistically significant flux is:
(2.26)
Where Nmin = minimum required number of aerosol to be observed over the measurement period, vd = the exchange velocity - see later section (mm s-1) and m = the time period of the measurement (s). The associated error in the calculated deposition velocity can be estimated from:
(2.27)
Where w = standard deviation in vertical wind velocity (m s-1) and Ntot = total number of aerosol observed during the measurement period. It can be shown, using the definition of deposition velocity (2.30) and (2.27), and assuming no extra error in the mean concentration, that the error in the flux is:
(2.28)
With rf = aerosol sample flow rate (cm3 s-1).
2.4.8 Eddy Covariance Limitations
There are several assumptions and requirements behind the eddy covariance technique. These are briefly reviewed here. Care must be taken in designing eddy covariance experiments, in analysing the data and in interpretation of the results with these factors in mind. The first of these is simply that in order to measure a flux attributable to a surface interaction, measurements must be made in the surface (constant flux) layer. This is something that must be addressed in the planning stage of an experiment.
2.4.9 Stationarity
It
is a fundamental assumption of the eddy covariance method that the
mean values of the quantities being measured do not change
significantly with time. Because of the Reynolds decomposition used
in the calculation, a changing concentration results in erroneous
values of the perturbations (x’, y’) being
used. The effect is shown in figure 2 for the covariance between a
steadily increasing and a stationary quantity. A subset of
instationarity errors are known as advection errors. These also
result from instationary values but are caused by shorter period
changes in concentration than those mentioned previously. Examples
include aerosol concentration as a polluted plume passes a
measurement location in a changing wind direction or a short duration
discharge within the footprint.
Instationarity as defined here can be corrected for by using linear de–trending (some workers use polynomial de–trending). There is no simple means of correcting for advection errors and they must accordingly be ‘planned out’ of experiments in advance as far as possible. As part of the normal quality control procedures used in this work, all power and co–spectra for aerosol measurements were inspected for behaviour significantly deviating from that expected (described later) and flagged as suspect as required.
2.4.10 Fetch and Homogeneity
Measurements made using the eddy covariance method are representative of a large area of the surface. This area is known as the fetch. In order for a measurement to be representative of a particular surface type or process, the whole of the fetch must have the same surface type. Since in theory the fetch extends an infinite distance upwind it is reasonable to arbitrarily select a proportion of the measured flux which should have come from the surface type of interest, usually around 90%. The need for the whole of the fetch to have the same surface type is known as the requirement for homogeneity.
Analytical solutions of the diffusion equation exist so that the relative contribution to the flux of various points upwind of the measurement point can be estimated. For example, Schuepp et. al. (1990) give:
(2.29)
Where Q = flux at the measurement point due to emission from one square metre of upwind surface (m-2 s-1), Q0 = a constant used to normalise the flux (m-2 s-1) and x = upwind distance (m). The Schuepp model is used in this work to estimate the contribution of certain fetch types to measured fluxes.
2.5 Interpretation of Fluxes
2.5.1 Deposition Velocity
When a species is depositing to a surface, its rate of deposition can be thought of as being controlled by the concentration of that species and by the rate at which it mixes downwards. The effective rate of downward movement of a species is called the deposition velocity, defined by:
(2.30)
Where vdm = measured deposition
velocity (conventionally mm s-1),
= average species concentration (e.g. cm-3 for
aerosol number) and f
= flux of trace species (e.g. for aerosol cm-2s-1).
C is a constant used to reconcile the units in the numerator
and denominator and to produce the desired units for vdm,
and the minus sign makes the convention that toward–surface
transport results in a positive deposition velocity. In fact, the
deposition velocity is higher than this measured vdm,
as gravitational settling is not measured by the eddy covariance
method. However, for 100 nm aerosol at STP, the gravitational
settling velocity (vg) predicted by the Stokes
equation is of the order of 1 m
s-1, so again the error in neglecting this has been
allowed in order to simplify the analysis. Therefore:
(2.31)
Where vd = actual deposition velocity at the measurement point (mm s-1). Deposition velocity is a useful tool because it is “semi–mechanistic”. It represents the effective velocity at which a trace substance moves towards a given surface type. It can then be used in conjunction with estimates of atmospheric concentration to infer deposition fluxes in modelling work. The interpretation of vdm is discussed further both here and in chapter 4 in the context of aerosol deposition modelling.
Analogous to the deposition velocity, but not widely used, is the emission velocity (ve). This is defined in exactly the same way as deposition velocity (2.30) but without the minus sign, so that transport away from the surface is positive by convention. ve is not as easily interpretable as vd, but can be useful. For example, it plays a part in the derivation of the aerosol concentration model described in chapter 6.
One tenable interpretation for emission velocity
in an emission situation comes from analogy with vd.
In (2.30), vd is the result of the process
f being
driven by
.
This is the way in which vd is used in some model
implementations. For emission velocity (in emission situations)
is the result of the process ve being driven by f.
This interpretation is discussed further and confirmed in chapter 6.
This being the case, ve is not useful in estimating
fluxes, but can be used as a crude empirical measure of vertical
ventilation of emitted species away from their source.
2.5.2 Resistance Analogy
Micrometeorological measurements of exchange of trace gases and aerosol are of little use unless they can be interpreted in a way that allows the effect of different surfaces to be accounted for as a way of generalising the results to other conditions. Another way of looking at deposition to surfaces is by analogy with electrical resistance. Taking Ohm’s law and replacing potential difference with concentration gradient (assuming zero concentration at the surface), and replacing current with flux gives:
(2.32)
Where rt = total resistance to deposition for species (s m-1). Comparing (2.30 – 2.32) it is clear that the resistance is equal to the inverse of deposition velocity for a given species and surface type (but note the difference in units). Often rt is split into two components to separate the turbulent transport to the top of the laminar sub–layer from the surface interaction. These are known as ra (aerodynamic resistance) and rs (surface resistance) respectively, such that (using the usual inversion rule for summation of series resistances):
(2.33)
In fact, in equation (2.33) the gravitational settling velocity is not accounted for, but as described previously this is negligibly small. The aerodynamic resistance term ra can be estimated, allowing the surface resistance to be calculated from (2.33). From Monteith and Unsworth (1990):
(2.34)
Which are both directly measured by the eddy covariance technique. The estimate of rs provided by (2.33 – 2.34) allows assessment of the surface capture efficiency, represented by the surface deposition velocity vds, where:
(2.35)
The surface deposition velocity or resistance can be further split into terms representing resistance to transport across the laminar boundary layer and resistance to deposition on different surface element types:
(2.36)
Where rb = resistance to transport across the laminar boundary layer or sub–layer and rc = canopy resistance, the resistance to deposition on a surface. Values of rb are not often directly determined in deposition studies (Monteith and Unsworth, 1990), but for “normal” (in crop micrometeorology terms) surfaces, the following has been proposed:
(2.37)
This approximation is thought to be adequate for 0.1 < u* < 0.5 m s-1 to calculate rb for heat and water vapour transfer to crops. Values from (2.37) would obviously need to be changed for species with molecular diffusivities significantly different to those effective for heat and water vapour transfer. An effective approximation for rb allows assessment of the canopy resistance rc to be made from deposition measurements. Reliable measurements of rc can, in some cases, allow identification of the component of a canopy responsible for uptake of a pollutant. For example, a gas which deposits rapidly during the day but not at night, and which appears to follow the same deposition trends as CO2 might be supposed to be depositing to plant stomata.
rc can be further subdivided into resistances to deposition on different canopy and surface elements. rg is the resistance of the ground surface, usually soil in canopy studies. The importance of rg obviously varies according to the density of the canopy and factors such as surface wetness for some species. rp denotes the resistance to surface uptake by chemical reaction, which is important for species such as ozone in cases where it reacts with chemicals on leaf surfaces. rs is the resistance to deposition on leaf surfaces (usually stomata, but a cuticular resistance term is sometimes added in parallel with rs) without chemical reaction. rs varies with time as plant stomata open and close according to ambient temperature, humidity and PAR (photosynthetically active radiation) availability. While in itself a useful quantity, rs can be divided further into resistances to deposition for different parts of the canopy surface. For example stomatal resistance and cuticular resistance are in general very different for most species. The amount and type of data required to perform an analysis in this depth is not available in field studies conducted here. Figure 3 shows a schematic of the different resistances, in series and parallel, which make up the total resistance to deposition rt.
Other techniques can assist with the interpretation of eddy covariance fluxes. Those used in this work include Fourier spectral and co–spectral analyses and probability density analysis. These are discussed briefly in the final sections of this chapter.
2.6 Fourier Analysis
2.6.1 Power Spectrum
Almost any data series can be adequately represented as a sum of sine and cosine waves of different frequencies and intensities. Fourier transforms allow the frequencies contributing to a series to be identified and the distribution of variance across those frequencies to be estimated. In this work only the forward Fourier transform is used to convert time series into frequency space, although an inverse transform does exist to convert from frequency space back to time series. The forward transform is defined, using Euler’s formula (cos(x) + sin(ix) = eix), with:
(2.35)
Where A(k) is a time series
consisting of N points (k = 0
N-1), nm = base period of the time series and n
= the series of frequencies equal to nm
k. The Fourier series is generally composed of a series of N
complex numbers. In micrometeorological applications the transform is
known as a power spectrum and is displayed as the modulus of the
complex Fourier series. In fact only half the series SA(n)
is displayed since the Fourier series is always symmetrical about the
point
and the second (redundant) half of the power series is always
removed.
fNy is called the Nyquist frequency and is very important in time series analysis. The reflection of the Fourier series about fNy is not a quirk in the definition of the transform. Rather it is a fundamental limit to the frequency information that can be extracted from a time series of finite frequency (length). If a wave of one–second period is to be resolved in a time series, the logging frequency of the measurement must be equal to at least 2 Hz (unless the instrument has a longer response time than the wave period, in which case it simply cannot be done). This limit manifests itself in frequency domain analysis as aliasing of signals faster than fNy into the transform at lower frequencies. This can directly affect flux measurements if an instrument is being sampled too infrequently. Given that the Nyquist frequency is exactly half the sampling frequency, a general guide can be derived for the minimum logging frequency for eddy covariance measurements using (2.3):
(2.36)
This means that the logging frequency used must be significantly faster than twice the dominant frequencies of motion predicted by mixing length theory as an absolute minimum.
2.6.2 Properties of the Atmospheric Power Spectrum
The properties of power spectra can be studied and
compared given the correct normalisations for the frequency and
variance axes. For comparison, power spectra are normalised by their
variances so that they integrate to unity and conventionally are
weighted by frequency. Frequency weighting gives the spectra the
units of variance (e.g. m2 s-2) rather
than variance per logarithmic frequency interval. This weighting
retains the feature of the log–log representation of the
unweighted spectra that power relationships appear as straight lines,
a distinct advantage in understanding the spectral dynamics of
turbulence. The frequency used for atmospheric applications is
non–dimensionalised by
to allow comparison of spectra taken at different heights and wind
speeds. For instrumentation testing (see chapters 3 and 4) the
natural frequency n (Hz) is used. For atmospheric applications
then, a typical power spectrum looks something like figure 3.
The two envelopes shown in figure 3 were suggested by Kaimal et. al. (1972) after an analysis of data from the Kansas AFCRL experiment (1968). It appears that under stable conditions the peak in the spectral intensity moves smoothly from lower to higher frequencies with increasing stability. Under unstable conditions they found that spectra tended to fall into the upper envelope, but that the behaviour of the peak was less ordered as a function of stability. This may indicate that in unstable conditions the peak frequency scales only with z (Kaimal and Finnigan, 1994). It may also be a result of the problems with surface layer scaling discussed in the section on Townsend’s Hypothesis.
The low frequency end of the spectrum is determined by large–scale motions, – the spectral range where energy enters the system. The +1 gradient observed at the lowest frequencies up to the spectral peak indicates that variance is low at low frequencies, and hence that energy is entering the cascade at scales up to the peak frequency. For the Kansas data, Kaimal et. al. (1972) estimate that the peak frequencies are approximately related by:
(2.37)
Where fm = peak in the frequency weighted power spectrum of species . The higher frequency part of the transform is very much simpler than this. Since no energy enters the cascade at high frequencies, the power spectrum is solely the result of movement of energy through vortex stretching and breaking processes to lower frequencies. Kolmogorov (1941) produced a dimensional analysis in which it was shown that the power spectrum in this region follows an f--5/3 power law (f-2/3 for the frequency–weighted spectrum). The magnitude of the inertial sub–range is governed by the power and frequency of the peak in the energy containing range. Kaimal and Finnigan (1994) suggest the following approximations:
(2.38)
(2.39)
(2.40)
Where is a function representing the dissipation rate of turbulent kinetic energy and h represents the variability of temperature. T* is a surface layer scaling temperature. They are defined by:
(2.41a)
(2.41b)
(2.42a)
(2.42b)
(2.43)
The gradients of these approximations for inertial sub–range spectra are used to assess the performance of the flux measurements during the GRAMINAE integrated experiment. At the highest frequency end of the inertial sub–range a power relation (gradient on the log–log axes used) of +1 is occasionally observed. This is not normally a real atmospheric phenomenon, but rather indicates that at that frequency range the time series is heavily influenced by white noise. This can be a good indication of the frequency response of an instrument, or can indicate that the fluctuations in that quantity are not resolvable by the instrument at that frequency range (figure 3). There are more sophisticated models for power spectra available than the ones shown above for the inertial sub–range, but these are not used in this thesis.
2.6.3 Co-spectral analysis
The co-spectrum is a broadly similar concept to the power spectrum. Rather than the variance in one time series at each frequency, the co-spectrum shows the covariance between two time series as a function of frequency. Usually co-spectra are calculated for the vertical wind speed perturbation (w’) with the value of a time series which has been used to calculate a flux by eddy covariance. This provides an estimate of the scales of motion (frequencies) contributing to a measured flux. The co-spectrum of two series A(k) and B(k), CAB(n) is defined as:
(2.44)
Where N is the sum of the number of samples in
A(k) and B(k) and the asterisk represents
the complex conjugate. Co-spectra are, similarly to power spectra,
often presented as
with non–dimensionalised frequency
.
Note the use of covariance (rAB) rather than
variance in the normalisation. This means that all such co-spectra
should have similar (order of magnitude) values, and spectra from
different heights and periods can be easily compared.
Frequency weighted co-spectra are presented in a manner analogous to the power spectrum representation presented in the previous section and shown in figure 3. There are also similar models available for co-spectra, and these will be introduced later as required. For now, it will simply be noted that the co-spectrum has a similar inertial sub–range to the power spectrum, and that this shows up as a f–4/3 power law in the frequency weighted plots.
2.7 Probability density analysis
Probability density analysis is used to display the joint two dimensional probability density function for turbulent fluctuations in two quantities. Data from one eddy covariance calculation period are plotted on each chart. Each axis represents the perturbation of one variable normalised by its standard deviation in order to allow intercomparison of different plots. Data are assigned to one of (in this work) eighty bins on the axes according to the perturbation value of each of the two series, making a total of 1600 bins. The frequency of data falling into each bin is displayed as a colour on the chart, showing the probability density multiplied by the number of samples used.
Probability density analysis is used to help identify the nature of a measured flux. For example it can be used to discriminate between the effects of strong vertical gradients in the measured property and strong vertical mixing. The plot shown in figure 4 is the result of a negative concentration gradient giving rise to net deposition. There are some apparent instantaneous emissions, but on average deposition is the dominant process. The analysis is used in the section on urban micrometeorology to illustrate the effects of instability on the aerosol flux. Figure 4 shows a hypothetical probability density and explains how it may be interpreted.