An analytic function for AOD

The volcanic radiative “forcing” produced by major stratospheric eruptions is usually estimated from observations of atmospheric optical depth ( AOD ). This is obtained by a variety of methods including ground and satellite based measurements. Several groups have worked on producing global averages or geographically gridded data. One such dataset is provided by NASA and covers the Mt Pinatubo eruption period.

Data is given as a series of layers of the stratosphere ( above 15km ) and does not cover the relatively minor contributions from the troposphere.

It is sometimes convenient to have a mathematical description of the time evolution of such a “forcing” rather than a time series of discrete data points.

It was found that the convolution of a exponential function with another exponential will fit Mt Pinatubo disturbance very closely.

The above graph selects the tropical bands from the data and sums the optical densities over the different height layers. Since AOD is a logarithmic measurement based on the transmission of light, summing different layer values is appropriate to get the total AOD for the column.

[ Tropical bands used: 19.6 11.7 3.9 -3.9 -11.7 -19.6 ]

The two time constants found to give a close fit were 3 and 9 months. Non-linear least squares fitting produced 2.95 and 9.1 months. It was found that adding an offset of 0.3 months from the time of eruptions gave a better fit. This can probably be explained as being due to the initial geographic dispersal of the ejecta after the eruption.

This curve provides convenient analytical mathematical description of the time evolution of AOD but it would be good to have an explanation that this does actually represent the physical processes and is not just a fortuitous resemblance of form.

The third line on the graph is a further convolution which represents a linear relaxation response of the climate system to the radiative perturbation. The derivation of the 8 month climate response is discussed here:

The major atmospheric effects of volcanoes, once the initial dust and ash has settled, is in the creation of aerosols. Sulphur dioxide combines with water and water vapour to produce a fine aerosol of dilute sulphuric acid. These persist in the stratosphere for several years.

In kinetic chemistry these are referred to as linear rate processes, where the rate of the process is determined by the concentration of reactant [1].

For a fixed initial amount, the concentration of the reactant decays exponentially with time. In engineering terms this can be called the impulse response of the system.

One way to calculate the system reaction to a varying input in such a situation, is by convolution of the continuous input with exponential impulse response.

Approximating the explosive injection of SO2 into the atmosphere as an instantaneous impulsion producing a step change in the amount of SO2 and applying a convolution to represent its conversion by reaction with water vapour into the dilute acid aerosol, then feeding this into a second convolution to represent the removal process, will produce a simple model of the evolution of aerosol concentration over time.

The analytical function resulting from these convolutions, where the reaction rates are different, is given by [2]:

λ1λ2 / ( λ1 – λ2 ) . ( exp1t – exp2t ) ; Eqn. 1

where λ1 and λ2 are the reciprocals of the time constants of the two reactions.

A special case where the two reactions have the same rate gives [3]:

λ2t . exp t ; Eqn. 2

Douglass & Knox 2005 [4] used a function, such as eqn. 2, to model AOD, though no explanation was given for the choice. It seems to have been chosen as a convenient function which provided a reasonably close fit to the data, rather than having a specific theoretical origin. Here is figure 2 from that paper with the authors’ original caption:

The single time-constant of this model was derived to 7.6 months, between the two values found above for the dual exponential model. It can be noted that this does not rise quickly enough, peaks a little too late and decays somewhat too quickly in the tail.

Though the fit is good and mean squared residuals are small, it does mean that two of key features that are of importance when applying this to studying the climate reaction are more notably affected: the overall scaling is less, and the timing of the peak is later.

Also, if the fitted AOD model decays too quickly, the implied climate response will be lengthened. This was another parameter that the paper sought to estimate.

None of these differences will be large and do not undermine the results of the paper, however the double decay model may provide a more accurate result and has the appeal of being grounded in a physical explanation of the chemical processes.

Investigation of the magnitude of the differences may justify using the single exponential model as a more parsimonious substitution for the physically grounded description in some situations. The impact on any derived parameters, as discussed above, and knock on effects to subsequent regression analyses and attribution studies should also be determined. A physically meaningful model would seem preferable in trying to understand the behaviour of the system, rather than just describe it in the most parsimonious way.


[1] U.C. Davis

[2] European Journal of Statistics and Probability
Vol.1, No.2, pp.1-8, December 2013
Oguntunde P.E1; Odetunmibi O.A2 ;and Adejumo, A. O 3.

[3] Dartmouth university course notes: section 7.2

[4] Douglass and Knox 2005: “Climate forcing by the volcanic eruption of Mount Pinatubo”