## Abstract

The explosion at the Buncefield oil depot in Hertfordshire, UK on 11 December 2005 produced the largest fire in Europe since the Second World War. The magnitude of the fire and the scale of the resulting plume thus present a stringent test of any mathematical model of buoyant plumes. A large-eddy simulation of the Boussinesq equations with suitable initial conditions is shown to reproduce the characteristics of the observed plume; both the height of the plume above the source and the direction of the downwind spread agree with the observations. This supports the use of the Boussinesq assumption, even for such a powerful plume as the one generated by the Buncefield fire. The presence of a realistic water vapour profile does not lead to significant additional latent heating of the plume, but does lead to a small increase in the final rise height of the plume due to the increased buoyancy provided by the water vapour. Our simulations include the interaction of radiation with the aerosol in the plume, and reproduce the observed optical thickness of the plume and the reduction of solar radiation reaching the ground. Far downwind of the source, solar radiation is shown to play a role in lofting the laterally spreading plume, but the manner in which it does so depends on the aerosol concentration. In the case of high aerosol concentration, the thickness of the plume increases; the incoming solar radiation is absorbed over such a small depth that only the top of the plume is lofted upwards and the level of maximum concentration remains almost unchanged relative to the case with no radiation. When the aerosol concentration is low, the whole plume is heated by the incoming solar radiation and the lofting is more coherent, with the result that the level of maximum concentration increases relative to the case with no radiation, but the thickness of the plume increases only slightly.

## 1. Introduction

The explosion at the Buncefield oil depot in Hertfordshire, UK on 11 December 2005 produced the largest peacetime fire in Europe to date. The explosion occurred at 06.00 and the resulting fire burnt for 4 days, with the bulk of the fuel burnt during the first 2 days. The incident was triggered by the ignition of approximately 300 tonnes of refined oil that had spilled over from a storage tank. Initial estimates of the fuel burnt over the 4-day event ranged from 7.1×10^{7} l (Mather *et al*. 2007) to 1.05×10^{8} l (Webster *et al*. 2006), which amounted to a release of thermal energy ranging from 2.5×10^{15} to 3.6×10^{15} J. Not surprisingly, there was considerable uncertainty in these estimates and later estimates of the amount of fuel burnt were more conservative, being of the order of 3.5–4.0×10^{7} l (Buncefield Investigation Board 2006). Owing to the efforts to contain and extinguish the fire, the rate at which the fuel was burnt decreased significantly after the second day. Observations of the resulting plume showed that it rose to approximately 3 km in height and then spread laterally in a fan-like shape to the southeast at lower levels and to the southwest at upper levels. Fortunately, the prevailing stable meteorological conditions ensured that very little of the emitted material reached the ground in a concentrated form. Observations of the plume showed that it consisted mostly of black carbon (soot) with unexceptional levels of other pollutants such as SO_{2}, CO, O_{3}, NO_{x} and hydrocarbons. The large amounts of black carbon in the plume suggest that the burning of the fuel was not very efficient and led, in part, to the lower estimates quoted above for the amount of fuel burnt. More detailed observations of the plume and its associated chemistry are described in Webster *et al*. (2006), Mather *et al*. (2007) and the reports by the Buncefield Investigation Board (2006).

The mathematical modelling of buoyant plumes has its origins in the physics of thermal convection. In the formulation of such models, it is commonplace to invoke the Boussinesq assumption, which states that if the temperature-induced change in the density of a fluid is sufficiently small, then, to a good approximation, any change in the density varies linearly with the temperature change. The result is that variations in density are only important when considering the buoyancy of the fluid and thus allows the fluid to be treated as incompressible. The aim of this paper is to assess the ability of large-eddy simulation (LES) of the Boussinesq equations to reproduce the characteristics of the plume generated by the Buncefield fire. LES is a commonly used approach to modelling and studying turbulent flows, in which the largest scales are resolved explicitly and the effects of the smallest scales are modelled. The advantage of this approach over a direct numerical simulation (DNS) of the governing equations, in which all scales are resolved, is its ability to simulate flows with very large Reynolds numbers, which are typical of the atmosphere and are beyond the capability of current DNSs, at relatively modest computational cost. In the context of atmospheric flows, LES of buoyant plumes has been considered by Nieuwstadt & de Valk (1987) and McGrattan *et al*. (1996). Computational resources place a strong constraint on domain size. A nested approach of an Eulerian chemistry transport model incorporating a simplified oil fire plume model has been used by Vautard *et al*. (2007) to simulate the dispersion of the Buncefield plume to hundreds of kilometres from the source. Large-scale dispersion of the Buncefield plume has also been performed using the Met Office Lagrangian stochastic model NAME (Webster *et al*. 2006, in press). Since our primary interest is in the near-source properties of the plume, such as its maximum height, a relatively small domain is sufficient for our purposes. Additionally, we incorporate the effects of moisture and radiation with a view to assessing their importance on the structure and evolution of the Buncefield plume.

The extremely large amount of energy released by the fire would almost certainly have resulted in large density differences and vertical accelerations of the order of the gravitational acceleration, which would violate the assumptions that underlie the Boussinesq approximation. However, Rooney & Linden (1997) and Woods (1997) have shown that these effects are only important in a region close to the source. For a point source, they define a length scale, *z*_{B}, above which the plume can be adequately described within the Boussinesq framework. It is given bywhere *F*_{B} is the buoyancy flux at the source; *g* is the acceleration due to gravity; and the entrainment constant, *α*, is taken to be 0.1 (this is the standard value of *α* and is consistent with the present model; see Devenish & Rooney (submitted) for details). For a buoyancy flux of 5.6×10^{5} m^{4} s^{−3}, commensurate with the energy release quoted above, this gives *z*_{B}≈300 m. Of course, a realistic plume has a finite source size and hence, for fixed *F*_{B}, the plume is likely to become more Boussinesq as the area of the source increases. Thus, *z*_{B} represents an upper bound on the region of non-Boussinesq effects. Since our primary interest is in the height attained by the plume and the subsequent lateral spread of material, the relatively small region over which the Boussinesq assumption is invalid does not concern us here. Moreover, LES has inherent difficulties in resolving the largest eddies in the vicinity of a surface. Finally, observations of the plume showed that its maximum height was of the order of 3 km, well below the scale height over which changes in the ambient density (induced by pressure changes) make the Boussinesq assumption invalid.

## 2. The governing equations

The Boussinesq equations for an incompressible fluid form the basis of the large-eddy model. The evolution of the resolved, or filtered, velocity, ** u**, is governed by(2.1)where

*p*is the pressure perturbation about the hydrostatic pressure;

*ρ*is a constant reference density; and

**is the unit vector in the vertical direction. The buoyancy is given by , where**

*k**θ*is the instantaneous potential temperature at time

*t*and position

**=(**

*x**x*,

*y*,

*z*) and

*θ*

_{0}is a reference potential temperature (here, taken to be 275 K, which is the observed temperature at the ground). The evolution of

*θ*is given by(2.2)The residual stress tensor,

**, and the heat flux,**

*σ***, represent the influence of the subgrid-scale (unresolved) terms. The Smagorinsky–Lilly model is used to approximate the subgrid terms in which the residual stress tensor takes the form , where and**

*H**i*,

*j*=1, …, 3. The subgrid eddy viscosity is given by , where

*λ*is the mixing length; ; and

*f*(

*x*) is a function of the Richardson number, , used to describe the stability of atmospheric flows. The form of

*f*reflects the observation that turbulence persists for

*Ri*≲1 and subsides for

*Ri*≳1. The mixing length is given by , where

*z*

_{0}is the roughness length;

*κ*=0.4 is the von Kármán constant; and (as in the standard Smagorinsky model), where

*Δ*is the grid spacing and

*C*

_{s}=0.23 is the Smagorinsky coefficient. For further details and the analogous model for the subgrid scalar flux, see Brown (1999) and references therein.

The momentum equations are solved using a centred-difference method based on Piacsek & Williams (1970). The evolution equation for the potential temperature (additional scalar fields are introduced in §4) is solved using a flux conservation method (Leonard *et al*. 1996). The Poisson equation for the pressure is solved using a fast Fourier transform. Periodic boundary conditions are imposed on the lateral boundaries for the velocity and pressure, but not for the scalar fields (the temperature perturbation, the scalar concentration and the moisture fields to be introduced in §4). On the lateral boundaries, the scalar fields are set to their initial values, except for the scalar concentration, which is set to zero. A sufficiently large domain is chosen so that the periodicity has little effect on the plume over the duration of the simulation. No-slip and free-slip boundary conditions are imposed, respectively, on the lower and upper boundaries of the domain. The surface stress is determined using Monin–Obukhov similarity theory and the top of the domain is stress free.

Initial conditions are chosen to be as consistent as possible with the observed meteorological conditions on the day of the explosion. The initial potential temperature profile, shown in figure 1*a* as a function of height, is taken from the Met Office operational general circulation model, the unified model (MetUM), and reflects the stable conditions of 11 December 2005. Comparison of temperature profiles over the course of the day did not show significant changes. The nearest radiosonde observations were made at Nottingham (approx. 120 km to the northwest) and Herstmonceux (approx. 100 km to the southeast) and, although they show some differences, are also indicative of stable conditions. For a more detailed discussion of the meteorological conditions on and after the day of the explosion, see Webster *et al*. (2007).

The initial velocity profile is also taken from the MetUM and is shown in figure 1*b*. Although the wind varies more than the temperature during the course of the day, as we will see, the lateral dispersion of the plume closely follows the wind profile, and so any changes to the wind profile would result in concomitant changes to the lateral dispersion. Since the synoptic situation did not vary significantly over the course of the day, we do not consider changes in the temperature and wind profiles any further. The wind profile in the LES does not change significantly from its initial state during the course of the simulation.

A surface heat flux is imposed in a circular source of radius 100 m in the centre of the lower boundary and is zero everywhere else. Based on the values quoted above for the energy released by the fire and the heating rates over the 4 days, an approximate estimate of the initial surface heat flux density is 500 000 W m^{−2} and the corresponding buoyancy flux density is 18 m^{2} s^{−3}. In calculating the buoyancy flux density, we have assumed standard values of the reference temperature and the specific heat capacity of air at constant pressure. According to Kaye & Laby (1995, p. 76), the specific heat capacity at constant pressure varies by 8.3 per cent between 275 and 775 K, which is sufficiently small to justify the assumption of constant specific heat capacity and is in the spirit of the Boussinesq assumption. A scalar flux is introduced in the plume source in order to aid visualization of the plume. In §5, the scalar field is identified with black carbon in order to study the effect of radiation on the plume. In this case, observations suggest an emission rate of 10–20 kg s^{−1} (Webster *et al*. 2006).

Unless stated otherwise, the horizontal domain is 10 ×10 km, with a resolution of 50 m in each direction. The vertical extent of the domain is 5 km, with a resolution of 50 m. Doubling either the horizontal or vertical resolution did not produce significantly different results; for a more detailed discussion of domain size and resolution, see Devenish & Rooney (submitted) and Devenish *et al*. (submitted). Where time averages are presented, these are based on statistics that are sampled every minute for 1 hour (after a short spin-up period). In the following, we are primarily interested in the height of the plume above the source and its subsequent lateral spread downwind of the source. In order to determine the former and gain a better understanding of the rising plume, we present statistics along the vertical centreline of the plume. As the horizontal mean wind will cause the centreline of the plume to meander around the geometric centre of the domain, wherever we present centreline statistics, these are calculated by joining the points on each horizontal plane at which the quantity concerned, for example the mean vertical velocity, is maximum.

## 3. Dry simulations

Figure 2*a* shows a typical instantaneous three-dimensional snapshot of the scalar field for an emission rate of 20 kg s^{−1}. The structure of the plume is clearly visible with an approximately conical form in the lower part and a laterally spreading plume in the upper region. The vertical extent of the plume and the fan-like spread of material from the plume in a southeast to south westerly direction with rising height can also be discerned from the image. The direction of the spread of material is consistent with the wind profile shown in figure 1*b*. Another feature to note is that the peak of the plume is higher than the laterally spreading material. The visual assessment of the image depends, of course, on which isosurfaces are chosen, and the range of concentrations on a given surface can be seen in the contour plots of the time-averaged concentration in figure 2*b–e*. Figure 2*b* is a vertical slice taken at 275 m east of the source, which shows the southward spread of the detrained material and that the maximum height of the plume is above the level at which most of the material is detrained. Figure 2*c–e* shows horizontal slices of the scalar field at 1500, 2000 and 2500 m and clearly demonstrates the effect of the wind turning with height. These results are broadly in agreement with the observations.

In a stably stratified environment, the density of a buoyant plume will increase with height due to entrainment, while that of the environment decreases with height. At some point, the density difference vanishes and, above this height, forces act on the plume to reduce its momentum and eventually bring it to rest. Assuming that the rate of entrainment is proportional to the vertical velocity, the amount of entrained fluid in this upper region is small and so the fluid in the topmost part of the plume will fall back some distance as it spreads out laterally (though not as far as the level at which the density difference first vanishes). This overshoot is visible in figure 2*a*,*b* and is also reflected in the different levels at which the vertical velocity and buoyancy (relative to the ambient potential temperature) become effectively negligible. The time-averaged centreline vertical velocity, , is shown in figure 3*a* and the time-averaged centreline buoyancy, , is shown in figure 3*b*. It is clear from figure 3*a*,*b* that the centreline velocity effectively vanishes significantly above the level at which the centreline buoyancy effectively vanishes and that both levels are consistent with the scalar field shown in figure 2*b*. Figure 3*a*,*b* shows that the root-mean-square (r.m.s.) centreline buoyancy is very close to the mean centreline buoyancy, whereas the r.m.s. centreline velocity increases relative to the mean near the top of the plume. In the region between the level of neutral buoyancy and the plume top, the mixing generated by the descending and ascending parts of the plume increases the amount of turbulence (relative to the lower part of the plume and to an equivalent plume in a uniform environment). For a more detailed discussion of buoyant plumes in a stably stratified environment including a comparison of LES with mathematical plume theory and sensitivity tests, see Devenish & Rooney (submitted) and Devenish *et al*. (submitted). The results of this section show that the mean field quantities computed from the LES data agree reasonably well with the observations and provide a benchmark against which the effects of moisture and radiation can be considered in §§4 and 5.

## 4. Simulations with moisture

A further series of simulations was conducted to assess the effect of a moist atmosphere on the development of the plume. Luderer *et al*. (2006) argued that latent heating due to condensation of entrained water vapour dominated the energy budget of the plume generated by the Chisholm forest fire in Canada in 2001. It is possible that latent heating may play a similar role in the Buncefield plume, although the good agreement between the results of §3 and the observations suggests that it will not be as important.

In this section, the governing equations of §2 are modified to take into account the effects of moisture on the flow, so that, in (2.1), *ρ* is now a dry reference density; ** u** is the velocity of dry air, which, to a good approximation, can be assumed to be equal to that of the moist air; and

*p*is the perturbation of the total (dry and moist) pressure about a dry hydrostatic pressure,

*p*

_{d}. The moist buoyancy is given by(4.1)where

*r*

_{v}is the water vapour mixing ratio;

*r*

_{l}is the liquid water mixing ratio; and

*ϵ*is the ratio of the molecular masses of dry air and water vapour. The thermodynamic equation (2.2) is replaced by(4.2)where D/D

*t*is the Lagrangian advective derivative;

*L*

_{v}is the latent heat of vaporization evaluated at a reference temperature

*T*

_{0}=273.15 K; and

*c*

_{p}is the specific heat capacity of dry air at constant pressure. In deriving the left-hand side of (4.2), we have assumed a reversible, adiabatic, process and that the total water mixing ratio is conserved, i.e. any liquid water stays within the cloud and there is no precipitation. The right-hand side then represents the diabatic heating. The subgrid model is modified to take into account the effects of moisture via a redefinition of the Richardson number (MacVean & Mason 1990). Two additional equations govern the evolution of the mixing ratios of vapour and liquid water,where is the subgrid scalar flux and is a source or sink of water vapour or liquid water. Liquid water condensate is produced whenever the water vapour mixing ratio is larger than the saturation mixing ratio, i.e. . Analytical expressions for

*r*

_{s}, which is a function of

*p*

_{d}and the local temperature,

*T*, can be derived from the Clausius–Clapeyron equation on making use of

*r*

_{s}=

*ϵe*

_{s}/

*p*

_{d}, where

*e*

_{s}is the saturation vapour pressure. If we assume that

*L*

_{v}depends linearly on temperature, then (Emanuel 1994, p. 116)(4.3)where

*c*

_{pv}is the specific heat capacity of water vapour at constant pressure;

*c*

_{l}is the specific heat capacity of liquid water; and

*R*

_{v}is the gas constant for water vapour. In deriving (4.3), we have neglected the temperature dependence of the specific heat capacities. For typical atmospheric temperatures, this assumption may be a reasonable first approximation, but for the temperatures in the Buncefield fire, it may no longer hold. However, very few data exist on how the specific heat capacities behave at very high temperatures and, as the pressure within the fire is unlikely to have reached sufficiently large values for liquid water to form, we ignore any dependence of the specific heat capacities on

*T*. Moreover, since the main effect of latent heating occurs near the top of the plume, where temperatures are much closer to typical atmospheric values, for our purposes, (4.3) remains appropriate.

The initial water vapour mixing ratio is taken from the MetUM and is shown in figure 1*c*. The initial liquid water mixing ratio is set to zero. Figure 3*c* shows the time-averaged centreline liquid water mixing ratio, . Not surprisingly, the maximum value of occurs near the top of the plume, where the temperature becomes cold enough for condensation to occur. The magnitude of is small compared with typical values for standard cloud types and suggests that latent heating due to condensation will be small. Figure 3*c* also shows considerable variation of about its time-averaged mean, indicating that water does not remain in its liquid state for long. We find that cloud only forms in the vicinity of, or within, the plume and is of limited spatial extent. Figure 3*a* compares the centreline velocity for the moist and dry simulations. The former is slightly larger than the latter, but the level at which both velocities effectively vanish is almost identical. This indicates that the additional heating due to the condensation of liquid water does not make a large contribution to the total heating of the plume. This is supported by the absence of any appreciable difference in the buoyancy (defined in terms of the potential temperature alone and relative to the ambient potential temperature) between the dry and moist simulations. A more important effect of moisture is the increase in the moist buoyancy, (4.1), which results from the water vapour. The plume entrains relatively moist air close to the ground and lifts it upwards so that at higher levels, where the ambient water vapour mixing ratio is lower (figure 1*c*), the plume has a higher concentration of water vapour relative to its surroundings. The result is that when the plume finally settles down to its level of neutral buoyancy, the additional buoyancy due to the water vapour means that this level is higher than for the corresponding dry plume. Figure 3*b* shows that the centreline moist buoyancy is slightly higher than the centreline buoyancy of the dry plume. In §5, we consider the variation of the mean height of the plume with downwind distance, in which the effect of the additional buoyancy due to the water vapour is more apparent.

The fire itself will produce water vapour as a by-product of combustion. Following Webster *et al*. (in press), if we assume that 3.5×10^{7} l of fuel were burnt, then this would amount to a release of 35 000 tonnes of water. If all this water vapour were to condense as liquid water, the resulting latent heat would be approximately 8.75×10^{13} J. This represents approximately 7 per cent of the total energy released by burning 3.5×10^{7} l of fuel. Thus, the addition of latent heating due to the condensation of water vapour released by the fire has only a small effect on the plume. It is estimated (Buncefield Investigation Board 2006) that 2.5×10^{7} l of water were added to the fire in an effort to contain and extinguish the fire. Assuming that all this water first evaporated and then recondensed near the top of the plume, an equivalent calculation shows that this would release a further 6.25×10^{13} J of energy, which is comparable with the value given above and so represents a relatively small percentage of the total heat released by the fire. In reality, of course, it is very unlikely that all the added water would evaporate and recondense, and so this estimate represents an upper bound.

## 5. Simulations with radiation

In this section, we consider the interaction of radiation with the plume. As argued by Herring & Hobbs (1994), radiation played a significant role in lofting the plumes generated by the Kuwaiti oil fires in 1991. Of course, these fires occurred in the summer and at a much lower latitude, and so there is no reason to expect radiative heating to play such a strong role for the Buncefield plume. However, it is of interest to assess the importance of radiation, particularly, far downwind of the source.

As the observations showed, the plume was formed mainly of black carbon, which is highly absorbent of radiation. Observations of the plume were made by the facility for airborne atmospheric measurements (FAAM) aircraft 48 hours after the explosion at approximately 78 km downwind of the source. Details of these measurements can be found in Webster *et al*. (2006, in press). Most of the aerosol within the plume had a typical diameter of 2.5 μm or less (PM_{2.5}) and the single scattering albedo was found to be 0.53 at 0.55 μm. This compares with 0.4 (at 0.55 μm) found by Mather *et al*. (2007) using ground-based measurements near the source. Since the single scattering albedo at a given wavelength, *λ*, is defined bywhere *k*_{s} is the mass scattering cross section (area scattered per unit aerosol mass) and *k*_{a} is the mass absorption cross section; *k*_{s} and *k*_{a} are approximately equal in this case. The mass extinction cross section, defined as , was also derived from the aircraft observations and hence allows *k*_{s} and *k*_{a} to be evaluated. Although *k*_{s} and *k*_{a} are approximately equal, because of multiple scattering, the plume as a whole is predominantly absorbing.

The radiative fluxes will be calculated using the two-stream equations, which we briefly review here. For a given wavelength, the optical depth for a layer of depth Δ*z* is given bywhere *ρ*_{a} is the aerosol concentration. The reflection coefficient for diffuse radiation, *R*_{λ}, is given bywhich follows from the solution of the two-stream equations using the practical improved flux method (Zdunkowski *et al*. 1980). Here, and , where *D* is the diffusivity factor. Note that, in the limit Δ*τ*_{λ}→∞, *R*_{λ} is independent of *D*. Furthermore, , where is the asymmetry of the aerosol. Since the scattering is predominantly in the forward direction, a typical value of is 3/4. Thus, in the limit of large optical depth, we get a value of *R*_{λ}≈0.04, which, despite the high opacity of the plume, means that very little of the incoming solar radiation is reflected. This is consistent with the observation that the value of *ω*_{λ} for the Buncefield plume is low compared with smoke from other sources such as biomass burning for which 0.8≲*ω*_{λ}≲0.95 (e.g. Mitchell *et al*. 2006; Johnson *et al*. 2008) and explains why the smoke from the Buncefield fire appeared black rather than grey or white.

As is commonplace in models of radiation in the atmosphere, short-wave (solar) and long-wave (terrestrial) radiation are treated separately. The short-wave region is divided into six wavelength bands and the long-wave region into eight bands. The solar irradiance at the top of the atmosphere, *S*_{0}, and the solar zenith angle are taken to be 1412.6 W m^{−2} and 74.97°, respectively (Liou 1980, p. 46), giving an insolation of 366.3 W m^{−2}. The long-wave radiation emitted by the surface is determined from the local temperature at the ground according to the Stefan–Boltzmann law. As well as interacting with the aerosol, the radiation also interacts with any liquid water, water vapour and other trace gases such as carbon dioxide and ozone. The heating, or cooling, rates associated with the absorption of radiation are determined from the radiative fluxes. These are calculated by summing the results of a number of quasi-monochromatic calculations in a vertical column. In the short-wave region, each quasi-monochromatic flux consists of an upward and downward diffuse flux and a direct solar flux. In the following, the downward and upward fluxes are both defined to be positive and the net flux is defined as the upward minus the downward flux. Mie scattering is used to calculate the scattered radiation by the aerosol particles, and Rayleigh scattering is used to calculate the scattering by air molecules. The effects of multiple scattering are also accounted for. The model surface reflects 10 per cent of the incoming short-wave radiation (consistent with observations), but there is no further interaction between the solar radiation and the surface and, consequently, no effect on the surface temperature or the surface heat flux. In the model, the terrestrial radiation also has no effect on the surface heat flux. More details of the radiation scheme can be found in Edwards & Slingo (1996) and Cusack *et al*. (1999).

We now consider the LES with the radiation scheme described above. Unless stated otherwise, all our simulations with radiation are for a plume in a moist atmosphere. Ideally, we would use an emission rate within the observed range (10–20 kg s^{−1}) and compare the LES results with the aircraft observations, but computational costs preclude a domain of 80 km (the distance from the source at which the aircraft observations were made) in each horizontal direction and a resolution of 50 m and so we cannot hope to reproduce the observations accurately. Instead, since the observed aerosol concentration at 78 km from the source was of the order of 5×10^{−8} kg m^{−3}, we scale the emission rate so that typical model aerosol concentrations near the downwind edge of the LES domain are of this order. We thus consider two emission rates: the maximum observed emission rate, *Q*_{1}=20 kg s^{−1}, and a scaled emission rate, *Q*_{2}=0.6 kg s^{−1}. In reality, of course, the aerosol concentration will reach the observed value in a plume that has been exposed to radiation for much longer and with a greater surface area than is possible in our simulations. To allow for this as much as possible, we use a domain of 10 km (*x*) by 20 km (*y*) with the same resolution. We specify a constant mean wind in the direction of the longer side of the domain whose magnitude is chosen so that the plume rises to approximately the same height as in the more realistic simulations. The source is placed on the midline of the domain 5 km from the upwind edge. In the rest of this section, we consider statistics of the aerosol concentration averaged in time and across the plume (in the *x* direction). In the following, we present the centre of mass of the plume, which is given bywhere *Χ*_{i} is the mean (in the above sense) aerosol concentration (mass per unit volume) in the *i*th layer; *z*_{i} is its height above the ground; and *N* is the number of vertical grid points. As the effects of radiation are most prominent in the region furthest from the source, wherever we present vertical profiles of the mean aerosol concentration, these are further averaged over the last 2500 m of the domain.

### (a) Winter simulations

Since the explosion and fire at the Buncefield oil depot occurred in winter, this will be the main focus of our investigations. Figure 4 shows the height of the centre of mass of the plume as a function of downwind distance. For reference, we also include the moist and dry plumes with no radiation and a *control case* in which only the water vapour and other trace gases interact with the radiation. Radiative heating has little impact on the maximum height of the plume above the source relative to the moist case with no radiation. More significant is the gradual increase in height with downwind distance, relative to the moist plume, of those plumes in which the aerosal interacts with the radiation relative to the moist plume. Note also that the centre of mass of the *Q*_{2} plume, , rises higher than that of the *Q*_{1} plume, . The corresponding vertical profiles of aerosol concentration are shown in figure 5. It shows that the maximum height attained by the *Q*_{1} plume is higher than that of the *Q*_{2} plume, but that the level of maximum concentration is higher for the *Q*_{2} plume than for the *Q*_{1} plume (which is almost unchanged relative to the moist plume).

The *Q*_{1} plume generates typical aerosol concentrations of approximately 10^{−6} kg m^{−3} near the downwind edge of the LES domain. This is broadly consistent with the near-source measurements of Mather *et al*. (2007). Although all of the incident solar radiation at the top of the laterally spreading plume, approximately 250 W m^{−2} at 2500 m, is absorbed by the model aerosol at this concentration, the radiative heating due to the short-wave radiation has little effect on the bulk of the plume and the level of maximum concentration. This apparently paradoxical situation can be understood as follows: taking a typical value of *k*_{a} in the short-wave region to be 3000 m^{2} kg^{−1}, the optical depth associated with an aerosol concentration of 10^{−6} kg m^{−3} is 0.15 for a layer of 50 m. The direct radiative flux density, *F*_{λ}, in the *n*th layer is given by , where *F*_{0} is the radiative flux density at the top of the plume and *μ* is the cosine of the solar zenith angle. Assuming that the aerosol concentration remains constant over a depth of 500 m, the radiative flux density decreases from 250 to 0.8 W m^{−2}; over 50 per cent of the incoming radiation is absorbed in the first 50 m alone. Indeed, the e-folding depth (the depth at which Δ*τ*_{λ}/*μ*=1) associated with these values of *k*_{a} and *ρ*_{a} is approximately 87 m. The radiative heating rates are given bywhere *ρ*_{d} is the density of dry air and Δ*F*_{λ} is the flux density divergence over a layer Δ*z*. The top layer has an average heating rate of approximately 8 K h^{−1}, which decreases to 0.04 K h^{−1} over the bottom (tenth) layer. Of course, the aerosol concentration varies with height, *k*_{e} varies with wavelength band and we have omitted the effects of scattering, and so these estimates cannot be considered precise, but they illustrate the rapid decrease in the heating rate from the top of the plume downwards. Since the lowest part of the plume is not warmed by the solar radiation, only the top-most part of the laterally spreading plume is lofted by the radiative heating.

For aerosol concentrations typical of the *Q*_{2} plume near the downwind edge of the LES domain, the optical depth associated with a 50 m layer is 7.5×10^{−3}. Assuming that this concentration remains constant across the depth of the plume, the penetration of the radiation is much greater in this case than for the *Q*_{1} plume above. Using the same value of *k*_{a} as above, the e-folding depth is approximately 1730 m (compared with 87 m above)! An equivalent calculation of the heating rates with Δ*z*=50 m over a depth of 500 m gives 0.5 K h^{−1} in the top layer and 0.4 K h^{−1} in the bottom layer. As above, these estimates cannot be considered precise, but they illustrate the greater uniformity of the heating rate over the depth of the plume when the aerosol concentration is low. Thus, although the heating rates at the top of the plume are smaller than those for high aerosol concentrations, all of the plume is heated by the incoming solar radiation and so the whole plume is lofted upwards. The actual heating rates as determined from the model will be discussed below. While this order of magnitude argument pertaining to the short-wave heating supports the observed trends in the LES data, namely that the bulk of the *Q*_{2} plume including the level of maximum concentration rises, that of the *Q*_{1} plume does not change appreciably, it does not explain why does not rise as high as . As we shall see in §5*c*, this is due to long-wave cooling.

We remark that the total optical aerosol thickness associated with the plume at the downwind edge of the LES domain is of the order of one for either of the two emission rates and is broadly consistent with the point measurements of Mather *et al*. (2007), both close to the source and 50 km downwind of the source. Furthermore, the aerosol mass loading measured by Mather *et al*. (2007) at 50 km downwind of the source lies between the aircraft measurements at 78 km from the source and the near-source measurements of Mather *et al*. (2007). Since our simulations are consistent with both the near-source measurements and the aircraft observations, we expect that a simulation with a sufficiently large domain to agree well with the measurements of Mather *et al*. (2007) at 50 km from the source.

Figure 6 shows the typical heating rates associated with short- and long-wave radiation for high (*Q*_{1}) and low (*Q*_{2}) aerosol concentrations near the downwind edge of the LES domain. Also shown for reference are the heating rates associated with the control case. Figure 6 shows that the heating rate associated with the absorption of short-wave radiation for high aerosol concentrations is not only larger than that for *Q*_{2} but also that its peak value occurs significantly higher. As shown in figure 7*b*, this is because the net (upward) short-wave flux density decreases much more rapidly with decreasing height for *Q*_{1} than for *Q*_{2}. Indeed, the net short-wave flux density at the ground is zero in the former case, whereas it is approximately 75 per cent of the control case when the emission rate is *Q*_{2}. (As the plume passed over an observation site at 50 km downwind of the source, Mather *et al*. (2007) measured a change in the clearness index, defined as the ratio of the total downward flux density to the solar irradiance at the top of the atmosphere, from 0.69 to 0.57. The clearness index for our control case is 0.7 and that for the *Q*_{2} plume is 0.52 (which corresponds to measurements made 78 km from the source). The larger decrease in our model results compared with that of Mather *et al*. (2007) closer to the source may be due to geometrical effects that are neglected in the one-dimensional radiation scheme used here.) As argued above, the rapid absorption of solar radiation when the aerosol concentration is high results in a large heating rate at the top of the *Q*_{1} plume, but negligible heating of the bulk of the plume, with the result that the level at which the concentration is largest remains almost unchanged relative to the moist plume. On the other hand, for the weakly absorbing case, the level at which the heating rate is maximum for the *Q*_{2} plume occurs at the level of maximum concentration of the moist plume and explains why the bulk of the *Q*_{2} plume is lofted upwards.

Since the typical absorption coefficient of aerosol for long-wave radiation is an order of magnitude smaller than that for short-wave radiation, the optical depth associated with long-wave radiation only becomes significantly large when the aerosol concentration is sufficiently high. The relatively strong long-wave cooling associated with *Q*_{1} is evident in figure 6, whereas that associated with *Q*_{2} does not differ significantly from the control case in which the aerosol is absent. The corresponding long-wave flux densities are shown in figure 7*a*. The maximum cooling occurs at the level of maximum concentration and is significantly larger than the control case until close to the top of the plume. The long-wave cooling associated with *Q*_{1} has a non-negligible impact on the net heating rate, particularly near the base of the plume, and, as we shall see, explains why does not rise as high as . In §5*c*, we consider the separate effects of long- and short-wave radiation on the plume in more detail.

### (b) Summer simulations

Although the explosion at the Buncefield oil depot took place in winter, the differences between the simulations with emission rates *Q*_{1} and *Q*_{2} are likely to be more pronounced for summertime conditions. (We ignore the possibility that summertime conditions could produce a well-mixed convective boundary layer that could significantly change the dynamics of the plume. However, if this boundary layer is shallow, the plume could penetrate the stable region above the capping inversion in which case the stable temperature profile we are assuming here may not be entirely unrealistic.) We conducted equivalent simulations to those described above in which only the solar zenith angle was changed to a value appropriate for the summer solstice, namely 28.33°. In contrast to the winter case, figure 4 shows that rises more rapidly with downwind distance than . However, figure 5*b* shows similar trends to those observed in figure 5*a*: there is more lofting of the top of the *Q*_{1} plume than for the equivalent winter case, but the level of maximum concentration is almost unchanged relative to the winter value. For the *Q*_{2} plume, the bulk of the plume including the level of maximum concentration has risen above the corresponding winter case. The top of the plume though has not risen significantly; an order of magnitude estimate of the heating rate in the top 50 m layer of this plume is 0.6 K h^{−1} (assuming that the incident radiation at the top of the plume is 1000 W m^{−2}), only 20 per cent higher than the corresponding winter estimate. As the solar elevation decreases, the length of the path traversed by the incident radiation increases. In the summer, *μ*≈0.9 and so, for a given optical depth, the decrease in the radiative flux density over a layer of 50 m is disproportionately small relative to the winter case. This compensating effect results from winter to summer ratios of *μ* and *F*_{0}, which are both close to 1 : 3.5. Indeed, since the heating rate is proportional to , to leading order, the heating rate is independent of the solar zenith angle.

### (c) Radiation sensitivities

In this section, we return to the winter simulations and consider the separate effects of long- and short-wave radiation, scattering and the role of water vapour and other trace gases. Figure 4 shows that the centre of mass of the control case introduced in §5*a* (where the radiation interacts with the moisture and other trace gases but not with the aerosol) does not vary significantly from that of the moist plume with no radiation. Thus, we may conclude that any change in height of the plume with downwind distance relative to the moist (or dry) plume is primarily due to the interaction of radiation with the aerosol. However, the complexity of the radiative transfer equations means that the combined radiative effects of both the aerosol and the trace gases (including water vapour) are not simply a linear superposition of the two. It is in this context that we investigate the effects of scattering on the heating rate. The additional scattering due to the aerosol could decrease the mean free path of the photons, with the result that a higher fraction of the short-wave radiation is absorbed by the water vapour and other trace gases. On the other hand, the increased level of scattering may simply result in more radiation being reflected to space. We conducted a number of tests to determine the effects of (Mie) scattering by the aerosol. For simplicity, we neglect the effects of Rayleigh scattering. As expected, figure 8*a* shows that the scattering is unimportant in the long-wave region for the *Q*_{1} plume. Relative to the cooling rates for the aerosol alone, the presence of water vapour and other trace gases reduces the cooling rate as the net loss of long-wave radiation from the plume to space is reduced. In the short-wave region (figure 8*b*) when no water vapour and other trace gases are present, the scattering by the aerosol reduces the heating rate because more radiation is scattered back to space. The addition of water vapour and other trace gases further reduces the heating rate, as some of the incident radiation is absorbed by these gases above the plume. These results indicate that, although the net effect of scattering by the aerosol is small, the scattering that does occur reduces the heating rate. The same trends are observed when the aerosol concentration is low, although the differences are much smaller.

We conducted a number of simulations to assess the relative importance of short- and long-wave radiation. The vertical profile of the aerosol concentration in the presence of *short-wave radiation only* is shown in figure 9. The results indicate that, as expected, the lofting of the plume is primarily a consequence of short-wave radiation and the differences between the high and low aerosol concentrations are consistent with those described above. The profiles for the simulations with both short- and long-wave radiation are shown for reference and indicate that long-wave cooling has a more significant impact when the aerosol concentration is high. Figure 10 shows the centres of mass of the plumes with and without long-wave radiation. The centres of mass of the two plumes with short-wave radiation only do not show significant differences between the high and low aerosol concentrations. In the *Q*_{1} plume, all the incoming short-wave radiation is absorbed, but only by the top of the plume where the aerosol concentration is relatively low. In the *Q*_{2} plume, on the other hand, some of the short-wave radiation penetrates down to the ground, but that which is absorbed occurs throughout the plume, including the region where the concentration is high. The result is that does not rise as high as may be expected, whereas rises higher than might be anticipated and in such a way that they are very similar. Also shown for reference are and with both short- and long-wave radiation. As may be anticipated from figure 9, differs more than between those plumes with and without long-wave radiation. This indicates that it is long-wave cooling that is primarily responsible for why does not rise as high as . It also explains why, for the summer simulations, is higher than (since the short-wave heating is stronger but the long-wave cooling remains almost unchanged relative to the winter simulations).

We now turn to the effects of long-wave radiation only. In the case of high aerosol concentration, the differences between the plumes with and without long-wave radiation suggest that long-wave cooling has an appreciable, if not dominant, effect. In an atmosphere with no water vapour or other trace gases, the radiation interacts only with the aerosol, and long-wave cooling reduces the buoyancy of the plume. This is clearly visible in the decrease (with downwind distance) of relative to the dry plume, as shown in figure 11*a*. The extent to which the buoyancy of the plume is reduced as a result of long-wave cooling increases with aerosol concentration. Thus, as shown in figure 11*b*, the centre of mass of the *Q*_{2} plume does not decrease much relative to that of the dry plume. Relative to the moist case, there is a slight decrease in the centre of mass of the *long-wave control plume* in which long-wave radiation interacts with the water vapour and other trace gases only (but *not* the aerosol). This is due to the cooling (relative to the ambient temperature) that results from an increased level of water vapour in the plume relative to its ambient value (as discussed in §4). The combined effect of the aerosol, water vapour and other trace gases with long-wave radiation shows a gradual decrease in the centre of mass relative to the moist case (and in the case of high aerosol concentration also the dry case), which is consistent with the cooling rate for a given aerosol concentration. Figure 12 shows that long-wave radiation has little impact on the vertical profile of aerosol concentration relative to the moist plume. There is a slight decrease in height of the *Q*_{1} plume, which is consistent with the decrease in the centre of mass relative to the moist plume over the last 2500 m of the LES domain.

It might be thought that long-wave cooling would generate entrainment by the plume of warm air above the plume, but the following scale analysis suggests that this is unlikely to play an important role. Assuming that buoyant fluid above the plume is mixed homogeneously throughout the plume, then the production of buoyancy within the plume is given by , where *w*_{e} is the entrainment velocity and *h* is the thickness of the plume. Assuming also that the plume is optically thin so that the cooling is uniform across the plume, then a resulting velocity scale is given by , where Δ*F* is the flux divergence across the plume. As we expect the creation and dissipation of buoyancy to be approximately in balance, we get and hence , where *α*_{R} is a constant of entrainment. Figure 7*a* suggests that the net long-wave flux divergence across the *Q*_{2} plume is approximately 50 W m^{−2} and, using a value of *α*_{R}=0.2 that is typically used for cloud-top entrainment in stratocumulus (e.g. Grenier & Bretherton 2001), we find that, for *h*=1000 m and a 5 K temperature rise over this height, *w*_{e}≈0.002 m s^{−1}. This suggests that ‘plume-top’ entrainment will not cause the plume to rise significantly over 1 hour. On the other hand, we note that between 2000 and 2500 m, where the top half of the plume is located, the vertical velocity variance is slightly larger for the *Q*_{2} plume than the *Q*_{1} plume, suggesting that there is more turbulence in the *Q*_{2} plume. However, it is worth emphasizing that the differences between the *Q*_{2} plume with long-wave radiation only and the moist plume are small (figure 12*b*), and that, in the broader context of the Buncefield plume, are merely subtleties. By far the most important quantity for modelling buoyant plumes is the surface heat flux, which determines the height of the plume above the source. Any subsequent lofting of the plume downwind of the source is primarily due to short-wave heating, and long-wave cooling plays only a secondary role.

## 6. Conclusions

We have shown that LES of the Boussinesq equations is able to simulate well the plume generated by the Buncefield fire. In particular, we were able to capture the maximum height of the plume, the level of neutral buoyancy and the fan-like structure of the laterally spreading material. This alone provides a justification for using the Boussinesq assumption in a situation where it might have been thought to be inappropriate and demonstrates that LES with a resolution of 50 m is capable of simulating the plume. However, it should be emphasized that the present model cannot describe the near-source region of such a strongly buoyant plume to a great degree of accuracy. With this caveat, the Boussinesq approximation provides an adequate description of the upper part of the plume, which is of most interest in making predictions of its dispersion. Of course, non-Boussinesq effects in the near-source region may have a cumulative effect on the upper part of the plume, but since the LES results agree broadly with the observations, we expect these effects to be small. The results presented here support the use of the Boussinesq assumption in dispersion models such as the Met Office Lagrangian stochastic model, NAME (Webster *et al*. 2006, in press), to model dispersion from large oil fires. A detailed comparison of the LES results with NAME can be found in Devenish *et al*. (submitted).

The LES results show that the additional heating due to the latent heat of condensation of water vapour has little effect on the maximum plume height above the source. A more important effect of moisture is the additional buoyancy provided by the water vapour with a consequent increase in the final level of neutral buoyancy relative to the dry plume. Radiation also has little effect on the maximum plume height above the source. However, far downwind of the source, the radiation acts to loft the laterally spreading plume, although the manner in which it does so depends on the aerosol concentration. In the case of high aerosol concentration, the thickness of the plume increases; the incoming solar radiation is absorbed over such a small depth that only the top of the plume is lofted upwards and the level of maximum concentration remains almost unchanged relative to the moist plume with no radiation. In the case of low aerosol concentration, the whole plume is heated by the incoming solar radiation and the lofting is more coherent, with the result that the level of maximum concentration increases relative to the moist plume, whereas the thickness of the plume increases only slightly. The solar radiation is primarily responsible for these effects, though long-wave cooling plays a small, counter-balancing, role when the aerosol concentration is high. We would expect that an equivalent simulation able to capture the spread of the plume to 80 km downwind of the source would first show lofting of the top of the plume (sufficiently far downwind of the source but within 15 km), but only further downwind would the level of maximum aerosol concentration begin to rise as a result of short-wave heating.

Finally, we speculate as to what would have happened if the Buncefield explosion had taken place at a different time or under different meteorological conditions. If the explosion had occurred at night, then, sufficiently far downwind of the source, some decrease in height of the plume might have been expected due to long-wave cooling. Summertime conditions would have increased the lofting of the plume due to radiative heating (as discussed in §5*b*). On the other hand, the convective conditions of a typical day in summer could have resulted in large amounts of emitted material being brought down to the ground. Rain would have led to increased levels of deposition of the emitted material. The presence of cloud could also have led to different radiative effects. Measurements by the FAAM aircraft showed that there was substantial cloud 2 days after the explosion and that the liquid water path associated with the cloud decreased significantly across the plume while the temperature of the plume increased. If our arguments in §5*c* about plume-top entrainment are correct, then the evaporation of the cloud would have been almost entirely due to short-wave heating.

## Acknowledgments

The authors would like to thank Nicolas Bellouin, James Haywood and Helen Webster for supplying data on the Buncefield fire and plume. © Crown Copyright 2008.

## Footnotes

- Received July 15, 2008.
- Accepted September 16, 2008.

- © 2008 The Royal Society