On Determination of Tropical Feedbacks

This article has been posted for discussion at Judith Curry’s “Climate Etc.”




Satellite data for the period surrounding the Mt Pinatubo eruption in 1991 provide a means of estimating the scale of the volcanic forcing in the tropics. A simple relaxation model is used to examine how temporal evolution of the climate response will differ from the that of the radiative forcing. Taking this difference into account is essential for proper detection and attribution of the relevant physical processes. Estimations are derived for both the forcing and time-constant of the climate response. These are found to support values derived in earlier studies which vary considerably from those of current climate models. The implications of these differences in inferring climate sensitivity are discussed. The study reveals the importance of secondary effects of major eruptions pointing to a persistent warming effect in the decade following the 1991 eruption.

The inadequacy of the traditional application of linear regression to climate data is highlighted, showing how this will typically lead to erroneous results and thus the likelihood of false attribution.

The implications of this false attribution to the post-2000 divergence between climate models and observations is discussed.

Keywords: climate; sensitivity; tropical feedback; global warming; climate change; Mount Pinatubo; volcanism; GCM; model divergence; stratosphere; ozone; multivariate; regression; ERBE; AOD ;


For the period surrounding the eruption of Mount Pinatubo in the Philippines in June 1991, there are detailed satellite data for both the top-of-atmosphere ( TOA ) radiation budget and atmospheric optical depth ( AOD ) that can be used to derive the changes in radiation entering the climate system during that period. Analysis of these radiation measurements allows an assessment of the system response to changes in the radiation budget. The Mt. Pinatubo eruption provides a particularly useful natural experiment, since the spread of the effects were centred in the tropics and dispersed fairly evenly between the hemispheres. [1] fig.6 ( This meta-study provides a lot of detail about the composition and dispersal of the aerosol cloud.)

The present study investigates the tropical climate response in the years following the eruption. It is found that a simple relaxation response, driven by volcanic radiative ‘forcing’, provides a close match with the variations in TOA energy budget. The derived scaling factor to convert AOD into radiative flux, supports earlier estimations [2] based on observations from the 1982 El Chichon eruption. These observationally derived values of the strength of the radiative disturbance caused by major stratospheric eruptions are considerably greater than those currently used as input parameters for general circulation climate models (GCMs). This has important implications for estimations of climate sensitivity to radiative ‘forcings’ and attribution of the various internal and external drivers of the climate system.


Changes in net TOA radiative flux, measured by satellite, were compared to volcanic forcing estimated from measurements of atmospheric optical depth. Regional optical depth data, with monthly resolution, are available in latitude bands for four height ranges between 15 and 35 km [DS1] and these values were averaged from 20S to 20N to provide a monthly mean time series for the tropics. Since optical depth is a logarithmic scale, the values for the four height bands were added at each geographic location. Lacis et al [2] suggest that aerosol radiative forcing can be approximated by a linear scaling factor of AOD over the range of values concerned. This is the approach usually adopted in IPCC reviewed assessments and is used here. As a result the vertical summations are averaged across the tropical latitude range for comparison with radiation data.

Tropical TOA net radiation flux is provided by Earth Radiation Budget Experiment ( ERBE ) [DS2].

One notable effect of the eruption of Mt Pinatubo on tropical energy balance is a variation in the nature of the annual cycle, as see by subtracting the pre-eruption, mean annual variation. As well as the annual cycle due to the eccentricity of the Earth’s orbit which peaks at perihelion around 4th January, the sun passes over the tropics twice per year and the mean annual cycle in the tropics shows two peaks: one in March the other in Aug/Sept, with a minimum in June and a lesser dip in January. Following the eruption, the residual variability also shows a six monthly cycle. The semi-annual variability changed and only returned to similar levels after the 1998 El Nino.

Figure 1 showing the variability of the annual cycle in net TOA radiation flux in tropics. (Click to enlarge)

It follows that using a single period as the basis for the anomaly leaves significant annual residuals. To minimise the residual, three different periods each a multiple of 12 months were used to remove the annual variations: pre-1991; 1992-1995; post 1995.

The three resultant anomaly series were combined, ensuring the difference in the means of each period were respected. The mean of the earlier, pre-eruption annual cycles was taken as the zero reference for the whole series. There is a clearly repetitive variation during the pre-eruption period that produces a significant downward trend starting 18 months before the Mt. Pinatubo event. Since it may be important not to confound this with the variation due to the volcanic aerosols, it was characterised by fitting a simple cosine function. This was subtracted from the fitting period. Though the degree to which this can be reasonably assumed to continue is speculative, it seems some account needs to be taken of this pre-existing variability. The effect this has on the result of the analysis is assessed.

The break of four months in the ERBE data at the end of 1993 was filled with the anomaly mean for the period to provide a continuous series.

Figure 2 showing ERBE tropical TOA flux adaptive anomaly. (Click to enlarge)

Figure 2b showing ERBE tropical TOA flux adaptive anomaly with pre-eruption cyclic variability subtracted. (Click to enlarge)

Since the TOA flux represents the net sum of all “forcings” and the climate response, the difference between the volcanic forcing and the anomaly in the energy budget can be interpreted as the climate response to the radiative perturbation caused by the volcanic aerosols. This involves some approximations. Firstly, since the data is restricted to the tropical regions, the vertical energy budget does not fully account for energy entering and leaving the region. There is a persistent flow of energy out of the tropics both via ocean currents and atmospheric circulation. Variations in wind driven ocean currents and atmospheric circulation may be part of any climate feedback reaction.

Secondly, taking the difference between TOA flux and the calculated aerosol forcing at the top of the troposphere to represent the top of troposphere energy budget assumes negligible energy is accumulated or lost internally to the upper atmosphere. Although there is noticeable change in stratospheric temperature as a result of the eruption, the heat capacity of the rarefied upper atmosphere means this is negligible in this context.

Figure 3 showing changes in lower stratosphere temperature due to volcanism. (Click to enlarge)

A detailed study on the atmospheric physics and radiative effects of stratospheric aerosols by Lacis, Hansen & Sato [2] suggested that radiative forcing at the tropopause can be estimated by multiplying optical depth by a factor of 30 W / m2.

This value provides a reasonably close match to the initial change in ERBE TOA flux. However, later studies[3], [4] , attempting to reconcile climate model output with the surface temperature record have reduced the estimated magnitude of the effect stratospheric aerosols. With the latter adjustments, the initial effect on net TOA flux is notably greater than the calculated forcing, which is problematic; especially since Lacis et al reported that the initial cooling may be masked by the warming effect of larger particles ( > 1µm ). Indeed, in order for the calculated aerosol forcing to be as large as the initial changes in TOA flux, without invoking negative feedbacks, it is necessary to use a scaling of around 40 W/m2. A comparison of these values is shown in figure 4.

What is significant is that from just a few months after the eruption, the disturbance in TOA flux is consistently less than the volcanic forcing. This is evidence of a strong negative feedback in the tropical climate system acting to counter the volcanic perturbation. Just over a year after the eruption, it has fully corrected the radiation imbalance despite the disturbance in AOD still being at about 50% of its peak value. The net TOA reaction then remains positive until the “super” El Nino of 1998. This is still the case with reduced forcing values of Hansen et al as can also be seen in figure 4.

Figure 4 comparing volcanic of net TOA flux to various estimations aerosol forcing. ( Click to enlarge )

The fact that the climate is dominated by negative feedbacks is not controversial since this is a pre-requisite for overall system stability. The main stabilising feedback is the Planck response ( about 3.3 W/m2/K at typical ambient temperatures ). Other feedbacks will increase or decrease the net feedback around this base-line value. Where IPCC reports refer to net feedbacks being positive or negative, it is relative to this value. The true net feedback will always be negative.

It is clear that the climate system takes some time to respond to initial atmospheric changes. It has been pointed out that to correctly compare changes in radiative input to surface temperatures some kind of lag-correlation analysis is required : Spencer & Braswell 2011[5], Lindzen & Choi 2011[6] Trenberth et al 2010 [7]. All three show that correlation peaks with the temperature anomaly leading the change in radiation by about three months.

Figure 5 showing climate feedback response to Mt Pinatubo eruption. Volcanic forcing per Lacis et al. ( Click to enlarge )

After a few months, negative feedbacks begin to have a notable impact and the TOA flux anomaly declines more rapidly than the reduction in AOD. It is quickly apparent that a simple, fixed temporal lag is not an appropriate way to compare the aerosol forcing to its effects on the climate system.

The simplest physical response of a system to a disturbance would be a linear relaxation model, or “regression to the equilibrium”, where for a deviation of a variable X from its equilibrium value, there is a restoring action that is proportional to the magnitude of that deviation. The more it is out of equilibrium the quicker its rate of return. This kind of model is common in climatology and is central of the concept of climate sensitivity to changes in various climate “forcings”.

dX/dt= -k*X ; where k is a constant of proportionality.

The solution of this equation for an infinitesimally short impulse disturbance is a decaying exponential. This is called the impulse response of the system. The response to any change in the input can found by its convolution with this impulse response. This can be calculated quite simply since it is effectively a weighted running average calculation. It can also be found by algebraic solution of the ordinary differential equation if the input can be described by an analytic function. This is the method that was adopted in Douglass & Knox 2005 [17] comparing AOD to lower tropospheric temperature ( TLT ).

The effect of this kind of system response is a time lag as well as a degree of low-pass filtering which reduces the peak and produces a change in the profile of the time series, compared to that of the input forcing. In this context linear regression of the output and input is not physically meaningful and will give a seriously erroneous value of the presumed linear relationship.

The speed of the response is characterised by a constant parameter in the exponential function, often referred to as the ‘time-constant’ of the reaction. Once the time-constant parameter has been determined, the time-series of the system response can be calculated from the time-series of the forcing.

Here, the variation of the tropical climate is compared with a linear relaxation response to the volcanic forcing. The magnitude and time-constant constitute two free parameters and are found to provide a good match between the model and data. This is not surprising since any deviation from equilibrium in surface temperature will produce a change in the long-wave Planck radiation to oppose it. The radiative Planck feedback is the major negative feedback that ensures the general stability of the Earth’s climate. While the Planck feedback is proportional to the fourth power of the absolute temperature, it can be approximated as linear for small changes around typical ambient temperatures of about 300 kelvin. This is effectively a “single slab” ocean model but this is sufficient since diffusion into deeper water below the thermocline is minimal on this time scale. This was discussed in Douglass & Knox’s reply to Robock[18]

It is this delayed response curve that needs to be compared to changes in surface temperature in a regression analysis. Regressing the temperature change against the change in radiation is not physically meaningful unless the system can be assumed to equilibrate much faster than the period of the transient being studied, ie. on a time scale of a month or less. This is clearly not the case, yet many studies have been published which do precisely this, or worse multivariate regression, which compounds the problem. Santer et al 2014 [8], Trenberth et al 2010 [7], Dessler 2010 b [9], Dessler 2011 [10] Curiously, Douglass & Knox [17] initially calculate the relaxation response to AOD forcing and appropriately regress this against TLT but later in the same paper regress AOD directly against TLT and thus find an AOD scaling factor in agreement with the more recent Hansen estimations. This apparent inconsistency in their method confirms the origin of the lower estimations of the volcanic forcing.

The need to account for the fully developed response can be seen in figure 6. The thermal inertia of the ocean mixed layer integrates the instantaneous volcanic forcing as well as the effects of any climate feedbacks. This results in a lower, broader and delayed time series. As shown above, in a situation dominated by the Planck and other radiative feedbacks, this can be simply modelled with an exponential convolution. There is a delay due to the thermal inertia of the ocean mixed layer but this is not a simple fixed time delay. The relaxation to equilibrium response introduces a frequency dependent phase delay that changes the profile of the time series. Simply shifting the volcanic forcing forward by about a year would line up the “bumps” but not match the profile of the two variables. Therefore neither simple regression nor a lagged regression will correctly associate the two variables: the differences in the temporal evolution of the two would lead to a lower correlation and hence a reduced regression result leading to incorrect scaling of the two quantities.

Santer et al 2014 attempts to remove ENSO and volcanic signals by a novel iterative regression technique. A detailed account, provided in the supplementary information[8b], reports a residual artefact of the volcanic signal.

The modelled and observed tropospheric temperature residuals after removal of ENSO and volcano signals, τ , are characterized by two small maxima. These maxima occur roughly 1-2 years after the peak cooling caused by El Chichon and Pinatubo cooling signals.

Figure XXX. Santer et al 2014 supplementary figure 3 ( panel D )
“ENSO and volcano signals removed”

This description matches the period starting in mid-1992, shown in figure 6 below, where the climate response is greater than the forcing. It peaks about 1.5 years after the peak in AOD, as described. Their supplementary fig.3 shows a very clear dip and later peak following Pinatubo. This corresponds to the difference between the forcing and the correctly calculated climate response shown in fig. 6. Similarly, the 1997-98 El Nino is clearly visible in the graph of observational data (not reproduced here) labeled “ENSO and volcano signals removed”. This failure to recognise the correct nature and timing of the volcanic signal leads to an incorrect regression analysis, incomplete removal and presumably incorrect scaling of the other regression variables in the iterative process. This is likely to lead to spurious attributions and incorrect conclusions.

Figure 6 showing tropical feedback as relaxation response to volcanic aerosol forcing ( pre-eruption cycle removed) ( Click to enlarge )

The delayed climatic response to radiative change corresponds to the negative quadrant in figure 3b of Spencer and Braswell (2011) [5] excerpted below, where temperature lags radiative change. It shows the peak temperature response lagging around 12 months behind the radiative change. The timing of this peak is in close agreement with TOA response in figure 6 above, despite SB11 being derived from later CERES ( Clouds and the Earth’s Radiant Energy System ) satellite data from the post-2000 period with negligible volcanism.

This emphasises that the value for the correlation in the SB11 graph will be under-estimated, as pointed out by the authors:

Diagnosis of feedback cannot easily be made in such situations, because the radiative forcing decorrelates the co-variations between temperature and radiative flux.

The present analysis attempts to address that problem by analysing the fully developed response.

Figure 7. Lagged-correlation plot of post-2000 CERES data from Spencer & Braswell 2011. (Negative lag: radiation change leads temperature change.)

The equation of the relationship of the climate response : ( TOA net flux anomaly – volcanic forcing ) being proportional to an exponential regression of AOD, is re-arranged to enable an empirical estimation of the scaling factor by linear regression.

TOA -VF * AOD = -VF * k * exp_AOD eqn. 1

-TOA = VF * ( AOD – k * exp_AOD ) eqn. 2

VF is the volcanic scaling factor to convert ( positive ) AOD into a radiation flux anomaly in W/m2. The exp_AOD term is the exponential convolution of the AOD data, a function of the time-constant tau, whose value is also to be estimated from the data. This exp_AOD quantity is multiplied by a constant of proportionality, k. Since TOA net flux is conventionally given as positive downwards, it is negated in equation 2 to give a positive VF comparable to the values given by Lacis, Hansen, etc.

Average pre-eruption TOA flux was taken as the zero for TOA anomaly and, since the pre-eruption AOD was also very small, no constant term was included.

Since the relaxation response effectively filters out ( integrates ) much of high frequency variability giving a less noisy series, this was taken as the independent variable for regression. This choice acts to minimise regression dilution due to the presence of measurement error and non-linear variability in the independent variable.

Regression dilution is an important and pervasive problem that is often overlooked in published work in climatology, notably in attempts to derive an estimation of climate sensitivity from temperature and radiation measurements and from climate model output. Santer et al 2014, Trenberth et al 2010, Dessler 2011, Dessler 2010b, Spencer & Braswell 2011.

The convention of choosing temperature as the independent variable will lead to spuriously high sensitivity estimations. This was briefly discussed in the appendix of Forster & Gregory 2006 [11] , though ignored in the conclusions of the paper.

It has been suggested that a technique based on total least squares regression or bisector least squares regression gives a better fit, when errors in the data are uncharacterized (Isobe et al. 1990). For example, for 1985–96 both of these methods suggest YNET of around 3.5 +/- 2.0 W m-2 K-1 (a 0.7–2.4 K equilibrium surface temperature increase for 2 ϫ CO2), and this should be compared to our 1.0–3.6 K range quoted in the conclusions of the paper.

Regression results were thus examined for residual correlation.


Taking the TOA flux, less the volcanic forcing, to represent the climatic reaction to the eruption, gives a response that peaks about twelve months after the eruption, when the stratospheric aerosol load is still at about 50% of its peak value. This implies a strong negative feedback is actively countering the volcanic disturbance.

This delay in the response, due to thermal inertia in the system, also produces an extended period during which the direct volcanic effects are falling and the climate reaction is thus greater than the forcing. This results in a recovery period, during which there is an excess of incoming radiation compared to the pre-eruption period which, to an as yet undetermined degree, recovers the energy deficit accumulated during the first year when the volcanic forcing was stronger than the developing feedback. This presumably also accounts for at least part of the post eruption maxima noted in the residuals of Santer et al 2014.

Thus if the lagged nature of the climate response is ignored and direct linear regression between climate variables and optical depth are conducted, the later extended period of warming may be spuriously attributed to some other factor. This represents a fundamental flaw in multivariate regression studies such as Foster & Rahmstorf 2013 [12] and Santer et al 2014 [8] , among others, that could lead to seriously erroneous conclusions about the relative contributions of the various regression variables.

For the case where the pre-eruption variation is assumed to continue to underlie the ensuing reaction to the volcanic forcing, the ratio of the relaxation response to the aerosol forcing is found to be 0.86 +/- 0.07%, with a time constant of 8 months. This corresponds to the value reported in Douglass & Knox 2005 [17] derived from AOD and lower troposphere temperature data. The scaling factor to convert AOD into a flux anomaly was found to be 33 W/m2 +/-11%. With these parameters, the centre line of the remaining 6 month variability ( shown by the gaussian filter ) fits very tightly to the relaxation model.

Figure 6 showing tropical feedback as relaxation response to volcanic aerosol forcing ( pre-eruption cycle removed) ( Click to enlarge )

If the downward trend in the pre-eruption data is ignored (ie its cause is assumed to stop at the instant of the eruption ) the result is very similar ( 0.85 +/-0.09 and 32.4 W/m2 +/- 9% ) but leads to a significantly longer time-constant of 16 months. In this case, the fitted response does not fit nearly as well, as can be seen by comparing figures 6 and 8. The response is over-damped: poorly matching the post-eruption change, indicating that the corresponding time-constant is too long.

Figure 8 showing tropical climate relaxation response to volcanic aerosol forcing, fitted while ignoring pre-eruption variability. ( Click to enlarge )

The analysis with the pre-eruption cycle subtracted provides a generally flat residual ( figure 9 ), showing that it accounts well for the longer term response to the radiative disruption caused by the eruption.

It is also noted that the truncated peak, resulting from substitution of the mean of the annual cycle to fill the break in the ERBE satellite data, lies very close to the zero residual line.

Figure 9 showing the residual of the fitted relaxation response from the satellite derived, top-of-troposphere disturbance. ( Click to enlarge )

Since the magnitude of the pre-eruption variability in TOA flux, while smaller, is of the same order as the volcanic forcing and its period similar to that of the duration of the atmospheric disturbance, the time-constant of the derived response is quite sensitive to whether this cycle is removed or not. However, it does not have much impact on the fitted estimation of the scaling factor ( VF ) required to convert AOD into a flux anomaly or the proportion of the exponentially lagged forcing that matches the TOA flux anomaly.

Assuming that whatever was causing this variability stopped at the moment of the eruption seems unreasonable but whether it was as cyclic as it appears to be, or how long that pattern would continue is speculative. However, approximating it as a simple oscillation seems to be more satisfactory than ignoring it.

In either case, there is a strong support here for values close to the original Lacis et al 1992 calculations of volcanic forcing that were derived from physics-based analysis of observational data, as opposed to later attempts to reconcile the output of general circulation models by re-adjusting physical parameters.

Beyond the initial climate reaction analysed so far, it is noted that the excess incoming flux does not fall to zero. To see this effect more clearly, the deviation of the flux from the chosen pre-eruption reference value is integrated over the full period of the data. The result is shown in figure 10.

Figure 10 showing the cumulative integral of climate response to Mt Pinatubo eruption. ( Click to enlarge )

Pre-eruption variability produces a cumulative sum initially varying about zero. Two months after the eruption, when it is almost exactly zero, there is a sudden change as the climate reacts to the drop in energy entering the troposphere. From this point onwards there is an ever increasing amount of additional energy accumulating in the tropical lower climate system. With the exception of a small drop, apparently in reaction to the 1998 ‘super’ El Nino, this tendency continues to the end of the data.

While the simple relaxation model seems to adequately explain the initial four years following the Mt Pinatubo event, this does not explain it settling to a higher level.


Concerning the more recent estimations of aerosol forcing, it should be noted that there is a strong commonality of authors in the papers cited here, so rather than being the work of conflicting groups, the more recent weightings reflect the result of a change of approach: from direct physical modelling of the aerosol forcing in the 1992 paper, to the later attempts to reconcile general circulation model (GCM) output by altering the input parameters.

From Hansen et al 2002 [4] ( Emphasis added. )


We illustrate the global response to these forcings for the SI2000 model with specified sea surface temperature and with a simple Q-flux ocean, thus helping to characterize the efficacy of each forcing. The model yields good agreement with observed global temperature change and heat storage in the ocean. This agreement does not yield an improved assessment of climate sensitivity or a confirmation of the net climate forcing because of possible compensations with opposite changes of these quantities. Nevertheless, the results imply that observed global temperature change during the past 50 years is primarily a response to radiative forcings.

Form section 2.2.2. Radiative forcing:

Even with the aerosol properties known, there is
uncertainty in their climate forcing
. Using our SI2000
climate model to calculate the adjusted forcing for a
globally uniform stratospheric aerosol layer with optical
depth t = 0.1 at wavelength l = 0.55 mm yields a forcing of
2.1 W/m2 , and thus we infer that for small optical depths

Fa (W/m2) ~ 21 tau

In our earlier 9-layer model stratospheric warming after El Chichon and Pinatubo was about half of observed values (Figure 5 of F-C), while the stratospheric warming in our current model exceeds observations, as shown below.

As the authors point out, it all depends heavily upon the assumptions made about size distribution of the aerosols used when interpreting the raw data. In fact the newer estimation is shown, in figure 5a of the paper, to be about twice the observed values following Pinatubo and El Chichon . It is unclear why this is any better than half observed values in their earlier work. Clearly the attributions are still highly uncertain and the declared uncertainty of +/-15% appears optimistic.

3.3. Model Sensitivity

The bottom line is that, although there has been some
narrowing of the range of climate sensitivities that emerge
from realistic models [Del Genio and Wolf, 2000], models
still can be made to yield a wide range of sensitivities by
altering model parameterizations.

If the volcanic aerosol forcing is underestimated, other model parameters will have to be adjusted to produce a higher sensitivity. It is likely that the massive overshoot in the model response of TLS is an indication of this situation. It would appear that a better estimation lies between these two extremes. Possibly around 25 W/m2.

The present study examines the largest and most rapid changes in radiative forcing in the period for which detailed satellite observations are available. The aim is to estimate the aerosol forcing and the timing of the tropical climate response. It is thus not encumbered by trying to optimise estimations of a range climate metrics over half a century by experimental adjustment of a multitude of “parameters”.

The result is in agreement with the earlier estimation of Lacis et al.

A clue to the continued excess over the pre-eruption conditions can be found in the temperature of the lower stratosphere, shown in figure 3. Here too, the initial disturbance seems to have stabilised by early 1995 but there is a definitive step change from pre-eruption conditions.

Noting the complementary nature of the effects of impurities in the stratosphere on TLS and the lower climate system, this drop in TLS may be expected to be accompanied by an increase in the amount of incoming radiation penetrating into the troposphere. This is in agreement with the cumulative integral shown in figure 10 and the southern hemisphere sea temperatures shown in figure 12.

NASA Earth Observatory report [13] that after Mt Pinatubo, there was a 5% to 8% drop in stratospheric ozone. Presumably a similar removal happened after El Chichon in 1982 which saw an almost identical reduction in TLS.

Whether this is, in fact, the cause or whether other radiation blocking aerosols were flushed out along with the volcanic emissions, the effect seems clear and consistent and quite specifically linked to the eruption event. This is witnessed in both the stratospheric and tropical tropospheric data. Neither effect is attributable to the steadily increasing GHG forcing which did not record a step change in September 1991.

This raises yet another possibility for false attribution in multivariate regression studies and in attempts to arbitrarily manipulate GCM input parameters to engineer a similarity with the recent surface temperature records.

With the fitted scaling factor showing the change in tropical TOA net flux matches 85% of the tropical AOD forcing, the remaining 15% must be dispersed elsewhere within the climate system. That means either storage in deeper waters and/or changes in the horizontal energy budget, ie. interaction with extra-tropical regions.

Since the model fits the data very closely, the residual 15% will have the same time dependency profile as the 85%, so these tropical/ex-tropical variations can also be seen as part of the climate response to the volcanic disturbance. ie. the excess in horizontal net flux, occurring beyond 12 months after the event, is also supporting restoration of the energy deficit in extra-tropical regions by exporting heat energy. Since the major ocean gyres bring cooler temperate waters into the eastern parts of tropics in both hemispheres and export warmer waters at their western extents, this is probably a major vector of this variation in heat transportation as are changes in atmospheric processes like Hadley convection.

Extra-tropical regions were previously found to be more sensitive to radiative imbalance that the tropics:


Thus the remaining 15% may simply be the more stable tropical climate acting as a buffer and exerting a thermally stabilising influence on extra-tropical regions.

After the particulate matter and aerosols have dropped out there is also a long-term depletion of stratospheric ozone ( 5 to 8% less after Pinatubo ) [13]. Thompson & Solomon (2008) [14] examined how lower stratosphere temperature correlated with changes in ozone concentration and found that in addition to the initial warming caused by volcanic aerosols, TLS showed a notable ozone related cooling that persisted until 2003. They note that this is correlation study and do not imply causation.

Figure 11. Part of fig.2 from Thompson & Solomon 2008 showing the relationship of ozone and TLS.

A more recent paper by Soloman [15] concluded roughly equal, low level aerosol forcing existed before Mt Pinatubo and again since 2000, similarly implying a small additional warming due to lower aerosols in the decade following the eruption.

Several independent data sets show that stratospheric aerosols have increased in abundance since 2000. Near-global satellite aerosol data imply a negative radiative forcing due to stratospheric aerosol changes over this period of about -0.1 watt per square meter, reducing the recent global warming that would otherwise have occurred. Observations from earlier periods are limited but suggest an additional negative radiative forcing of about -0.1 watt per square meter from 1960 to 1990.

The values for volcanic aerosol forcing derived here being in agreement with the physics-based assessments of Lacis et al. imply much stronger negative feedbacks must be in operation in the tropics than those resulting from the currently used model “parameterisations” and the much weaker AOD scaling factor.

These two results indicate that secondary effects of volcanism may have actually contributed to the late 20th century warming. This, along with the absence of any major eruptions since Mt Pinatubo, could go a long way to explaining the discrepancy between climate models and the relative stability of observational temperatures measurements since the turn of the century.

Once the nature of the signal has been recognised in the much less noisy stratospheric record, a similar variability can be found to exist in southern hemisphere sea surface temperatures. The slower rise in SST being accounted for by the much larger thermal inertia of the ocean mixed layer. Taking TLS as an indication of the end of the negative volcancic forcing and beginning of the additional warming forcing, the apparent relaxation to a new equilibrium takes 3 to 4 years. Regarding this approximately as the 95% settling of three times-constant intervals ( e-folding time ) would be consistent with a time constant of between 12 and 16 months for extra-tropical southern hemisphere oceans.

These figures are far shorter than values cited in Santer 2014 ranging from 30 to 40 months which are said to characterise the behaviour of high sensitivity models and correspond to typical IPCC values of climate sensitivity. It is equally noted from lag regression plots in S&B 2011 and Trenberth that climate models are far removed from observational data in terms of producing correct temporal relationships of radiative forcing and temperature.

Figure 12. Comparing SH sea surface temperatures to lower troposphere temperature.

Though the initial effects are dominated by the relaxation response to the aerosol forcing, both figures 9 and 12 show a additional climate reaction has a significant effect later. This appears to be linked to ozone concentration and/or a reduction in other atmospheric aerosols. While these changes are clearly triggered by the eruptions, they should not be considered part of the “feedback” in the sense of the relaxation response fitted here. These effects will act in the same sense as direct radiative feedback and shorten the time-constant by causing a faster recovery. However, the settling time of extra-tropical SST to the combined changes in forcing indicates a time constant ( and hence climate sensitivity ) well short of the figures produced by analysing climate model behaviour reported in Santer et al 2014 [8].

IPCC on Clouds and Aerosols:

IPCC AR5 WG1 Full Report Jan 2014 : Chapter 7 Clouds and Aerosols:

No robust mechanisms contribute negative feedback.

The responses of other cloud types, such as those associated with deep convection, are not well determined.

Satellite remote sensing suggests that aerosol-related invigoration of deep convective clouds may generate more extensive anvils that radiate at cooler temperatures, are optically thinner, and generate a positive contribution to ERFaci (Koren et al., 2010b). The global influence on ERFaci is

Emphasis added.

WG1 are arguing from a position of self-declared ignorance on this critical aspect of how the climate system reacts to changes in radiative forcing. It is unclear how they can declare confidence levels of 95%, based on such an admittedly poor level of understanding of the key physical processes.


Analysis of satellite radiation measurements allows an assessment of the system response to changes in radiative forcing. This provides an estimation of the aerosol forcing that is in agreement with the range of physics-based calculations presented by Lacis et al in 1992 and is thus brings into question the much lower values currently used in GCM simulations.

The considerably higher values of aerosol forcing found here and in Lacis et al imply the presence of notably stronger negative feedbacks in tropical climate and hence imply a much lower range of sensitivity to radiative forcing than those currently used in the models.

The significant lag and ensuing post-eruption recovery period underlines the inadequacy of simple linear regression and multivariate regression in assessing the magnitude of various climate ‘forcings’ and their respective climate sensitivities. Use of such methods will suffer from regression dilution, omitted variable bias and can lead to seriously erroneous attributions.

Both the TLS cooling and the energy budget analysis presented here, imply a lasting warming effect on surface temperatures triggered by the Mt Pinatubo event. Unless these secondary effects are recognised, their mechanisms understood and correctly modelled, there is a strong likelihood of this warming being spuriously attributed to some other cause such as AGW.

When attempting to tune model parameters to reproduce the late 20th century climate record, an incorrectly small scaling of volcanic forcing, leading to a spuriously high sensitivity, will need to be counter-balanced by some other variable. This is commonly a spuriously high greenhouse effect, amplified by presumed positive feedbacks from other climate variables, less well constrained by observation ( such as water vapour and cloud ). In the presence of the substantial internal variability, this can be made to roughly match the data while both forcings are present ( pre-2000 ). However, in the absence of significant volcanism there will be a steadily increasing divergence. The erroneous attribution problem, along with the absence of any major eruptions since Mt Pinatubo, could explain much of the discrepancy between climate models and the relative stability of observational temperature measurements since the turn of the century.


[1] Self et al 1995
“The Atmospheric Impact of the 1991 Mount Pinatubo Eruption”

[2] Lacis et al 1992 : “Climate Forcing by Stratospheric Aerosols

[3] Hansen et al 1997 “Forcing and Chaos in interannual to decadal climate change”

[4] Hansen et al 2002 : “Climate forcings in Goddard Institute for Space Studies SI2000 simulations”

[5] Spencer & Braswell 2011: “On the Misdiagnosis of Surface Temperature Feedbacks from Variations in Earth’s Radiant Energy Balance”

[6] Lindzen & Choi 2011: “On the Observational Determination of Climate Sensitivity and Its Implications”

[7] Trenberth et al 2010: “Relationships between tropical sea surface temperature and top‐of‐atmosphere radiation”

 [8] Santer et al 2014: “Volcanic contribution to decadal changes in tropospheric temperature”

 [8b] Supplementary Information:

 [9] Dessler 2010 b “A Determination of the Cloud Feedback from Climate Variations over the Past Decade”

 [10] Dessler 2011 “Cloud variations and the Earth’s energy budget”

 [11]Forster & Gregory 2006
“The Climate Sensitivity and Its Components Diagnosed from Earth Radiation Budget Data”

[12] Foster and Rahmstorf 2011: “Global temperature evolution 1979-2010″

[13] NASA Earth Observatory

 [14] Thompson & Soloman 2008: “Understanding Recent Stratospheric Climate Change”

 [15] Susan Soloman 2011: “The Persistently Variable “Background” Stratospheric Aerosol Layer and Global Climate Change”

 [16] Trenberth 2002: “Changes in Tropical Clouds and Radiation”

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

 [18] Douglass & Knox 2005b: “Reply to comment by A. Robock on ‘‘Climate forcing by the volcanic eruption of Mount Pinatubo’’”

Data sources:
Derived from SAGE (Stratospheric Aerosol and Gas Experiment) instrument.
“The estimated uncertainty of the optical depth in our multiple wavelength retrievals [Lacis et al., 2000] using SAGE observations is typically several percent. “[4]

[DS2] ERBE TOA data: Earth Radiation Budget Experiment

data description:


[*] Explanatory notes:
Negative feedback is an engineering term that refers to a reaction that opposes its cause. It is not a value judgement about whether it is good or bad. In fact, negative feedbacks are essential in keeping a system stable. Positive feedbacks lead to instability and are thus generally “bad”.

Convolution is a mathematical process that is used, amongst other things, to implement digital filters. A relaxation response can be implemented by convolution with an exponential function. This can be regarded as an asymmetric filter with a non-linear phase response. It produces a delay relative to the input and alters it’s shape.

Regression dilution refers to the reduction in the slope estimations produced by least squares regression when there is significant error or noise in the x-variable. Under negligible x-errors OLS regression can be shown to produce the best unbiased linear estimation of the slope. However, this essential condition is often over-looked or ignored, leading to erroneously low estimations of the relationship between two quantities. This is explained in more detail here:

Supplementary Information

Detailed method description.

The break of four months in the ERBE data at the end of 1993 was filled with the anomaly mean for the period to provide a continuous series. This truncates what would have probably been a small peak in the data, marginally lowering the local average, but since this was not the primary period of interest this short defect was considered acceptable.

The regression was performed for a range of time-constants from 1 to 24 months. The period for the regression was from just before the eruption, up to 1994.7, when AOD had subsided and the magnitude of the annual cycle was found to change, indicating the end of the initial climatic response ( as determined from the adaptive anomaly shown in figure 2 ).

Once the scaling factor was obtained by linear regression for each value of tau, the values were checked for regression dilution by examining the correlation of the residual of the fit with the AOD regressor while varying the value of VF in the vicinity of the fitted value. This was found to give a regular curve with a minimum correlation that was very close to the fitted value. It was concluded that the least squares regression results were an accurate estimation of the presumed linear relationship.

The low values of time-constant resulted in very high values of VF ( eg. 154 for tau=1 month ) that were physically unrealistic and way out side the range of credible values. This effectively swamps the TOA anomaly term and is not meaningful.

tau = 6 months gave VF = 54 and this was taken as the lower limit for the range time constant to be considered.

The two free parameters of the regression calculations will ensure an approximate fit of the two curves for each value of the time constant. The latter could then be varied to find the response that best fitted the initial rise after the eruption, the peak and the fall-off during 1993-1995.

This presented a problem since, even with adaptive anomaly approach, there remains a significant, roughly cyclic, six month sub-annual variability that approaches in magnitude the climate response of interest. This may be at least in part a result aliasing of the diurnal signal to a period close to six months in the monthly ERBE data as pointed out by Trenberth [16]. The relative magnitude of the two also varies depending on the time-constant used, making a simple correlation test unhelpful in determining the best correlation with respect to the inter-annual variability. For this reason it could not be included as a regression variable.

Also the initial perturbation and the response, rise very quickly from zero to maximum in about two months. This means that low-pass filtering to remove the 6 month signal will also attenuate the initial part of the effect under investigation. For this reason a light gaussian filter ( sigma = 2 months ) was used. This has 50% attenuation at 2.3 months. comparable to a 4 month running mean filter ( without the destructive distortions of the latter). This allowed a direct comparison of the rise, duration of the peak and rate of fall-off in the tail that is determined by the value of the time-constant parameter, and thus the value producing the best match to the climate response.

Due to the relatively course granularity of monthly data that limits the choice of the time-constant values, it was possible to determine a single value that produced the best match between the two variables.

The whole process was repeated with and without subtraction of the oscillatory, pre-eruption variability to determine the effects of this adjustment.

The uncertainty values given are those of the fitting process and reflect the variability of the data. They do not account for the experimental uncertainty in interpreting the satellite data.

Triple running mean filter

The following script will call a simple running mean three times with appropriate window size to do effect triple running mean, as described in the article (as amended for the asymmetric kernel to minimise negative leakage):


It requires the runmean.awk script found here:


Select code with mouse to copy elsewhere.


# call runmean.awk three times to compose triple running mean
# usage: ./r3m.sh file window_len   ; default window is 12 data point

if [ "x$1" == "x"  ]; then echo "$0 : err no file name,   usage: $0 filename "
 exit 1;
else fn=$1

if [ "x$2" == "x"  ]; then win=12; else win=$2; fi

#win2=`awk "BEGIN{print $win/ 1.3371}"`
#win3=`awk "BEGIN{print $win/ 1.3371/ 1.3371}"`

# asymmetric stages with following window ratios minimise neg. lobe in freq. response
k=1.15; k2=1.58

win2=`awk "BEGIN{print $win/ "$k"}"`
win3=`awk "BEGIN{print $win/ "$k2" }"`

outfile=`echo $fn | awk '{ print substr($1,1,length($1)-4) }'`

echo  "# triple running mean :  $win  $win2  $win3 " > $outfile

# echo $fn; echo $win; echo $win2; echo $win3

cat $fn | ./runmean.awk - $win |  ./runmean.awk - $win2 |  ./runmean.awk - $win3 >> $outfile

#echo "# cat $fn | ./runmean.awk - $win |  ./runmean.awk - $win2 |  ./runmean.awk - $win3 >> $outfile"

echo "# triple running mean :  $win  $win2  $win3 "
echo "# outfile = "$outfile;

On inappropriate use of OLS


Figure 1 showing conventional and inverse OLS linear regression fits on some real, observed climate variables.

Least squares regression is a very useful technique, widely used in almost all branches of science. It is usually one of the first techniques that is taught in schools for analysing experimental data.

It is also a technique that is misapplied almost as often as it is used correctly.

It can be shown that under certain conditions, the least square linear regression is the best estimation of the true relationship that can be derived from the available data. In statistics this is often called the ‘best, unbiased linear estimator’ of the slope. (Some, who enjoy contrived acronyms, abbreviate this to “BLUE”.)

However, it has been known since at least 1878 that this will under-estimate the slope if there is measurement error in the x-variable. (R. J. Adcock) [0]

There are two main conditions for this result to be an accurate estimation of the slope. One is that the deviations of the data from the true relationship are ‘normally’ or gaussian distributed. That is to say that they are of a random nature. This condition can be violated by significant periodic components in the data or excessive number of out-lying data points. The latter may often occur when only a small number of data points is available and the noise, even if random in nature, is not sufficiently sampled to average out.

The other main condition is that there be negligible error ( or non-linear varibility ) in the x variable. If this condition is not met, the OLS result derived from the data will almost always under-estimate the slope of the true relationship. This effect is sometimes referred to as regression dilution. The degree by which the slope is under-estimated is determined by the nature of the x and y errors but most strongly by those in x since they are required to be negligible for OLS to give the best estimation.

In this discussion “errors” can be understood to be both observational inaccuracies and any variability due to some factor other than the supposed linear relationship that it is sought to determine by linear regression of the two variables.

In certain circumstances regression dilution can be corrected for, but in order to do so, some knowledge of the nature and size of the both x and y errors has to be known. Typically this is not the case beyond knowing whether the x variable is a ‘controlled variable’ with negligible error, although several techniques have been developed to estimate the error in the estimation of the slope [7].

A controlled variable can usually be attained in a controlled experiment, or when studying a time series, provided that the date and time of observations have been recorded and documented in a precise and consistent manner. It is typically not the case when both sets of data are observations of different variables, as is the case when comparing two quantities in climatology.

One way to demonstrate the problem is to invert the x and y axes and repeat the OLS fit. If the result were valid, irrespective of orientation, the first slope would be the reciprocal of second one. However, this is only the case when there is very small errors in both variables, ie. the data is highly correlated ( grouped closely around a straight line ). In the case of one controlled variable and one error prone variable, the inverted result will be incorrect. In the case of two datasets containing observational error, both results will be wrong and the correct result will generally lie somewhere in between.

Another way to check the result is by examining the cross-correlation between the residual and the independent variable ( model – y vs x ), then repeat for incrementally larger values of the fitted ratio. Depending on the nature of the data, it will often be obvious that the OLS result does not represent the minimum correlation between the ordinate and the regressor, ie. it does not optimally account for co-variability of the two quantities.

In the latter situation, the two regression fits can be taken as bounding the likely true value but some knowledge of the relative errors is needed to decide where in that range the best estimation lies. There are a number of techniques such as bisecting the angle, taking the geometric mean (square root of the product), or some other average, but ultimately, they are no more objective unless driven by some knowledge of the relative errors. Clearly bisection would not be correct if one variable had low error, since the true slope would then be close to the OLS fit done with that quantity on the x axis.

Figure 2. A typical example of linear regression of two noisy variables produced from synthetic randomised data. ( Click to enlarge graph and access code to reproduce data and graph. )

Figure 2b. A typical example of correct application of linear regression to data with negligible x-errors. The regressed slope is very close to the true value, so close as to be indistinguishable visually. ( Click to enlarge )

The larger the x-errors, the greater the skew in the distribution and the greater the dilution effect.

An Illustration: the Spencer simple model.

The following case is used to illustrate the issue with ‘climate-like’ data. However, it should be emphasised that the problem is an objective mathematical one, the principal of which is independent of any particular test data used. Whether the following model is an accurate representation of climate ( it is not claimed to be ) has no bearing on the regression problem.

In a short article on his site Dr. Roy Spencer provided a simple, single-slab ocean, climate model with a predetermined feedback variable built into it. He observed that attempting to derive the climate sensitivity in the usual way consistently under-estimated the know feedback used to generate the data.

By specifying that sensitivity (with a total feedback parameter) in the model, one can see how an analysis of simulated satellite data will yield observations that routinely suggest a more sensitive climate system (lower feedback parameter) than was actually specified in the model run.

And if our climate system generates the illusion that it is sensitive, climate modelers will develop models that are also sensitive, and the more sensitive the climate model, the more global warming it will predict from adding greenhouse gasses to the atmosphere.

This is a very important observation. Regressing noisy radiative flux change against noisey temperature anomalies does consistently produce incorrectly high sensitivity. However, it is not an illusion created by the climate system, it is an illusion created by the incorrect application of OLS regression. Where there is error on both variables, the OLS slope is no longer an accurate estimation of the underlying linear relationship being sought.

Dr Spencer was kind enough to provide an implementation of the simple model in the form of a spread sheet download so that others may experiment and verify the effect.

To demonstrate this problem, the spreadsheet provided was modified to duplicate the dRad vs dTemp graph but with the axes inverted, ie. using exactly the same data for each run but additionally displaying it the other way around. Thus the ‘trend line’ provided by the spreadsheet is calculated with the variables inverted. No changes were made to the model.

Three values were used for the feedback variable were used in turn, two values: 0.9 and 1.9 that Roy Spencer suggests represent the range of IPCC values and 5.0 which he proposes as a value closer to that which he has derived from satellite observational data.

Here is a snap-shot of the spreadsheet showing a table of results from nine runs for each feedback parameter value. Both the conventional and the inverted regression slopes and their geometric mean have been tabulated.

Figure 3. Snap-shot of spreadsheet, click to enlarge.

Firstly this confirms Roy Spencer’s observation that the regression of dRad against dTemp consistently and significantly under-estimates the the feedback parameter used to create the data (and hence over-estimates climate sensitivity of the model). In this limited test, error is between a third and a half of the correct value. There is only one value of the conventional least squares slope that is greater than the respective feedback parameter value.

Secondly, it is noted that the geometric mean of the two OLS regressions does provide a reasonably close to the true feedback parameter, for the value derived from satellite observations. Variations are fairly evenly spread either side: the mean is only slightly higher than the true value and the standard deviation is about 9% of the mean.

However, for the two lower feedback values, representing the IPCC range of climate sensitivities, while the usual OLS regression is substantially less than the true value, the geometric mean over-estimates and does not provide a reliable correction over the range of feedbacks.

All the feedbacks represent a net negative feedback ( otherwise the climate system would be fundamentally unstable ). However, the IPCC range of values represents less negative feedbacks, thus a less stable climate. This can be seen reflected in the degree of variability in data plotted in the spreadsheet. The standard deviations of the slopes are also somewhat higher. This is not unexpected with less feedback controlling variations.

It can be concluded that the ratio of the proportional variability in the two quantities changes as a function of the degree of feedback in the system. The geometric mean of the two slopes does not provide a good estimation of the true feedback for the less stable configurations which have greater variability. This is in agreement with Isobe et al 1990 [7] which considers the merits of several regressions methods.

The simple model helps to see how this relates to Rad / Temp plots and climate sensitivity. However, the problem of regression dilution is a totally general mathematical result and can be reproduced from two series having a linear relationship with added random changes, as shown above.

What the papers say

A quick review of several recent papers on the problems of determining climate sensitivity shows a general lack of appreciation of the regression dilution problem.

Dessler 2010 b [1] :

Estimates of Earth’s climate sensitivity are uncertain, largely because of uncertainty in the long-term cloud feedback.

Spencer & Braswell 2011 [2] :

Abstract: The sensitivity of the climate system to an imposed radiative imbalance remains the largest source of uncertainty in projections of future anthropogenic climate change.

There seems to be agreement that this is the key problem in assessing future climate trends. However, many authors seem unaware of the regression problem and much published work on this issue seems to rely heavily on the false assumption that OLS regression of dRad against dTemp can be used to correctly determine this ratio, and hence various sensitivities and feedbacks.

Trenberth 2010 [3] :

To assess climate sensitivity from Earth radiation observations of limited duration and observed sea surface temperatures (SSTs) requires a closed and therefore global domain, equilibrium between the fields, and robust methods of dealing with noise. Noise arises from natural variability in the atmosphere and observational noise in precessing satellite observations.

Whether or not the results provide meaningful insight depends critically on assumptions, methods and the time scales ….

Indeed so. Unfortunately they then go on to contradict earlier work by Lindzen and Choi that did not rely on OLS regression by relying on inappropriate use of regression.

Spencer and Braswell 2011 [2]

As shown by SB10, the presence of any time-varying radiative forcing decorrelates the co-variations between radiative flux and temperature. Low correlations lead to regression-diagnosed feedback parameters biased toward zero, which corresponds to a borderline unstable climate system.

This is an important paper highlighting the need to use lagged regression to avoid the decorrelating effect of delays in the response. However, does not deal with the further attenuation due to regression dilution. It is ultimately still based on regression of two error laden-variables and thus does not recognise regression dilution that is also present in this situation. Thus it is likely that this paper is still over-estimating sensitivity.

Dessler 2011 [4] :

Using a more realistic value of σ(dF_ocean)/σ(dR_cloud) = 20, regression of TOA flux vs. dTs yields a slope that is within 0.4% of lamba.

Then in the conclusion of the paper, emphasis added:

warming). Rather, the evolution of the surface and atmosphere during ENSO variations are dominated by oceanic heat transport. This means in turn that regressions of TOA fluxes vs. δTs can be used to accurately estimate climate sensitivity or the magnitude of climate feedbacks.

Also from a previous paper:

Dessler 2010 b [1]

The impact of a spurious long-term trend in either dRall-sky or dRclear-sky is estimated by adding in a trend of T0.5 W/m 2/ decade into the CERES data. This changes the calculated feedback by T0.18 W/m2/K. Adding these errors in quadrature yields a total uncertainty of 0.74 and 0.77 W/m2/K in the calculations, using the ECMWF and MERRA reanalyses, respectively. Other sources of uncertainty are negligible.

The author is apparently unaware of the inaccuracy of regressing two uncontrolled variables is a major source of uncertainty and error.

Lindzen & Choi 2011 [5]

[Our] new method does moderately well in distinguishing positive from negative feedbacks and in quantifying negative feedbacks. In contrast, we show that simple regression methods used by several existing papers generally exaggerate positive feedbacks and even show positive feedbacks when actual feedbacks are negative.

… but we see clearly that the simple regression always under-estimates negative feedbacks and exaggerates positive feedbacks.

Here the authors have clearly noted that there is a problem with the regression based techniques and go into quite some detail in quantifying the problem, though they do not explicitly identify it as being due to the presence of uncertainty in the x-variable.

The L&C papers, to their credit, recognise that regression seriously under-estimates the slope and utilise other techniques to more correctly determine the ratio.

It seems the latter authors are exceptional in looking at the sensitivity question without relying on inappropriate use of linear regression. It is certainly part of the reason that their results are considerably lower than almost all other authors on this subject.

Forster & Gregory 2006 [8]

For less than perfectly correlated data, OLS regression of Q-N against δTs will tend to underestimate Y values and therefore overestimate the equilibrium climate sensitivity (see Isobe et al. 1990).

Another important reason for adopting our regression model was to reinforce the main conclusion of the paper: the suggestion of a relatively small equilibrium climate sensitivity. To show the robustness of this conclusion, we deliberately adopted the regression model that gave the highest climate sensitivity (smallest Y value). It has been suggested that a technique based on total least squares regression or bisector least squares regression gives a better fit, when errors in the data are uncharacterized (Isobe et al. 1990). For example, for 1985–96 both of these methods suggest YNET of around 3.5 +/- 2.0 W m2 K-1 ( a 0.7–2.4-K equilibrium surface temperature increase for 2 ϫ CO2 ), and this should be compared to our 1.0–3.6-K range quoted in the conclusions of the paper.

Here, the authors explicitly state the regression problem and its effect on the results of their study on sensitivity. However, when writing in 2005, they apparently feared that it would impede the acceptance of what was already a low value of climate sensitivity if they presented the mathematically more accurate, but lower figures.

It is interesting to note that Roy Spencer, in a non peer reviewed article, found an very similar figure of 3.66 W/m2/K by comparing ERBE data to MSU derived temperatures following Mt Pinatubo. [10]

So Forster and Gregory felt constrained to bury their best estimation of climate sensitivity, and the discussion of the regression problem in an appendix. In view of the ‘gatekeeper’ activities revealed in the Climategate emails, this may have been a wise judgement in 2005.

Now, eight years after the publication of F&G 2006, proper application of the best mathematical techniques available to correct this systematic over-estimation of climate sensitivity is long overdue.



Inappropriate use of linear regression can produce spurious and significantly low estimations of the true slope of a linear relationship if both variables have significant measurement error or other perturbing factors.

This is precisely the case when attempting to regress modelled or observed radiative flux against surface temperatures in order to estimate sensitivity of the climate system.

In the sense that this regression is conventionally done in climatology, it will under-estimate the net feedback factor (often denoted as ‘lambda’). Since climate sensitivity is defined as the reciprocal of this term, this results in an over-estimation of climate sensitivity.

If an incorrect evaluation of climate sensitivity from observations is used as a basis for the choice of parametrised inputs to climate models, the resulting model will be over sensitive and produce exaggerated warming. Similarly faulty analyses of their output will further inflate the apparent model sensitivity.

This situation may account for the difference between regression-based estimations of climate sensitivity and those produced by other methods. Many techniques to reduce this effect are available in the broader scientific literature, thought there is no single, generally applicable solution to the problem.

Those using linear regression to assess climate sensitivity need to account for this significant source of error when supplying uncertainly values in published estimations of climate sensitivity or take steps to address this issue.

A mathematical explantion of the origin of regression dilution is provided here:
On the origins of regression dilution



 [1] Dessler 2010 b “A Determination of the Cloud Feedback from Climate Variations over the Past Decade”


 [2] Spencer and Braswell 2011: “On the Misdiagnosis of Surface Temperature Feedbacks from Variations in Earth’s Radiant Energy Balance”


 [3] Trenberth et al 2010 “Relationships between tropical sea surface temperature and top‐of‐atmosphere radiation”


 [4] Dessler 2011
“Cloud variations and the Earth’s energy budget”


 [5] Lindzen & Choi 2001 “On the Observational Determination of Climate Sensitivity and Its Implications”


 [6] Nic Lewis : “A Sensitive Matter: How The IPCC Buried Evidence Showing Good News About Global Warming ”

 [7] Isobe et al 1990 “Linear Regression in Astronomy I”


Forster & Gregory 2006
“The Climate Sensitivity and Its Components Diagnosed from Earth Radiation Budget Data”


Adcock, R.J., 1878. A Problem in Least Squares. The Analyst, 5, 53-54.

Quirino Paris 2004 : “Robust Estimators of Errors-In-Variables Models Part I”


R Spencer “Revisiting the Pinatubo Eruption as a Test of Climate Sensitivity”


Gaussian low-pass filter

Click to enlarge graph.

The graph is scaled to unit frequency (eg 1 Hz or one cycle per year …).

It illustrates how a 159 day gaussian will filter daily data. For monthly data 12*159/365=5.23 , so the nearest would be sigma=5mo . However since that is slightly short of the intended value and gaussian leaks a little, since it is never fully zero, six months would be better to suppress a 12mo cycle.

Code ( use copy / paste within code block ).

#!/bin/awk -f

# pass input through 3 sigma gaussian filter where sigma, if not given, is 2 data points wide
# usage : ./gauss.awk filename <sigma=2> <scale_factor=1>
# optional scale_factor simply scales the output
# use  OFMT="%6.4f"
# sigma can be compared to the period of the -3dB point of the filter
# result is centred, ie not shift. dataset shortened by half window each end

# operates on files with two columns of numeric data
#  eg date as integer year or floating point decimal
# check whether data is continuous !!
# nov 2012  OFMT="%10.8f"

# FWHM (-3dB) = 2*sigma*sqrt(2*ln2)
# half height freq = sigma * 1.1774
# 3 sigma gives 99.7% accuracy
# eg. sigma=6 months can be used to remove most of an annual signal
#  though other filters with a zero pass may be a better choice

BEGIN { OFMT="%10.8f"
# ARGV[1]=filename; argv[0] is script name, hence ARGC>=1

  if (ARGV[1]=="") { print " # usage : ./gauss.awk filename <sigma=2> <scale_factor=1>" ; exit 1}

  pi= 3.14159265359811668006

  if  ( ARGC >3 ) {scaleby=ARGV[3];ARGV[3]=""} else {scaleby=1};
  if  ( ARGC >2 ) {sigma=ARGV[2];ARGV[2]=""} else {sigma=2};

  print "# filtering "ARGV[1]" with gaussian of sigma= ",sigma

  gw=3*sigma-1;  # gauss is approx zero at 3 sigma, use 3 sig window

# calculate normalised gaussian coeffs
  for (tot_wt=j=0;j<=gw;j++) {tot_wt+=gwt[-j]=gwt[j]=exp(-j*j/two_sig_sqr)/ root2pi_sigma};
  for (j=-gw;j<=gw;j++) {gwt[j]/=tot_wt};

# strip off last .xxx part of file name
# improve this  (doesn't work with paths like ../filename)

  if (ARGV[1]=="-") {
   else {
    print "# ",out_file >out_file;


( ($0 !~ /^#/) && ($0 !~ /^[" ",\t]*$/)  ) {

  if (ln>2*gw)
    for (j=-2*gw;j<=0;j++) {gauss+=ydata[ln+j]*gwt[j+gw]}
    if (out_file=="-")
      print xdata[ln-gw],gauss
    else {
      print NR,xdata[ln-gw],gauss
      print xdata[ln-gw],gauss >> out_file;
#    print  $1,$2;


  if (ARGV[1]=="") { exit 1}
  print "# gausssian window width = "gw+gw+1",done"
  print "# output file = "out_file

# for (i=0;i<=gw;i++){print gwt[i]}

Lanczos high-pass filter

The graph shows the frequency response of a 3-lobe Lanczos high-pass filter with central frequencies of 461 and 269 days, respectively.

These are chosen such that a periodicity of 1 year is either fully included or fully excluded in the filter’s pass-band. ie the 269 day filter will pass only frequencies shorter than 365 days and fully block annual cycles. 461 day filter will ensure that everything including the annual signal is unattenuated and excludes anything longer than about 1.6 years.

The equivalent low-pass filter is detailed here:


frequency axis is in units of “per year” not “year” as labelled.

Code: ( select text in code block and use copy/paste to get the code )

#!/bin/awk -f

#  usage: ./lanczos_HP.awk filename [-3dB_period=15] [lobes=3]"

# purpose: apply a high-pass filter by convolvution of 2nd column of file with lanczos windowed sinc function
# column 1 data should be evenly spaced without gaps.
# optional second parameter is period of -3dB attenuation, as a period expressed  in data intervals
#    default = 15 ; eg. 15 monthly periods : half power at 15mo, close to 100% pass at 12 mo.
#                   to exclude 12mo, use 9 . This will pass 6mo and higher freq unattenuated.
#                   for daily data use 461 and 269 days , respectively.
# optional third parameter is the number of lobes in impulse response. More gives sharper transition but more frequency distortion.
#    default=3; 3 and 4 are good.  ### >=5 may produce significant ringing artefacts ###

function sinc(x) {
  if (x==0) return 1; else return (sin(x)/x);

BEGIN{  OFMT = "%8.6f"


  if ( ARGC >3 ) {lobes=ARGV[3];ARGV[3]=""} else {lobes=5};
  if ( ARGC >2 ) {period=ARGV[2];ARGV[2]=""} else {period=15};
  if ( ARGC <2) { print "### usage: ./lanczos_HP.awk filename [-3dB_period=15] [lobes=3]"; err=1; exit 1; }
  else {len=ARGV[2]; ARGC=2;}

  pi= 3.14159265359811668006 ;
  kwid= 2*halfk;         # kernel = -halfk....+halfk = kwid+1 pts ie odd number
  w_taper=twopi/kwid;    # first lanczos lobe -pi to +pi to fill window
  w=twopi/period;        # sinc passes zero pi, 2*pi, 3*pi ....

# calculate normalised kernel coeffs
  for (tot_wt=j=0;j<=halfk;j++) {tot_wt+=gwt[-j]=gwt[j]=sinc(j*w)*sinc(j*w_taper)};
  for (j=-halfk;j<=halfk;j++) {gwt[j]/=tot_wt};

# for HP filter, invert LP and add unity impulse:
  for (j=-gw;j<=gw;j++) {gwt[j]/=-tot_wt}

# strip off last .xxx part of file name
# improve this  (doesn't work with paths like ../filename)

  if (ARGV[1]=="-") {
  else {
    print "# ",out_file >out_file;
    print "# ",out_file

  ln=-1;  # init linenumber counter

($0 !~ /^#/)&&($0 != "") {    # strip out comments, headers and blank lines

  if (ln>2*halfk)
    for (j=-2*halfk;j<=0;j++) {sum+=ydata[ln+j]*gwt[j+halfk]}
    if (out_file=="-")
      print xdata[ln-halfk],sum
    else {
      print NR,xdata[ln-halfk],sum
      print xdata[ln-halfk],sum >> out_file;
  } #endif ln


    if (err) exit err;
    if (FNR<=kwid) print " ### insufficient data to fill "kwid+1" point filter buffer"; exit 1;

    print "# "lobes"-lobe lanczos kernel width = "2*halfk+1",done"

    print "# output file = "out_file
# uncomment following to dump kernel (impluse response):
# sum=0; for (j=-halfk;j<=halfk;j++) {sum+= gwt[j]; printf ("%.4f\n", gwt[j])}
# print "# abs sum of coeffs = ",sum;


Lanczos low-pass filter

Click to enlarge graph.

frequency axis is in units of “per year” not “year” as labelled.

Equivalent high-pass filter:


Discussion: https://climategrog.wordpress.com/?attachment_id=659

Source Code: ( to get code, use browser copy / paste function within code block below).

#!/bin/awk -f

#  usage: ./lanczos_LP.awk filename [50%_attn_period=18] [lobes=3]"

# purpose: convolve 2nd column of file with lanczos windowed sinc function
# this provides low-pass filter with a fast transition from pass-band to stop-band , with minimal ripple in both.
# column 1 data should be evenly spaced without gaps.
# optional second parameter is period of 50% attenuation, as a period expressed in data intervals
# optional third parameter is the number of lobes in impulse response.
# zero pass req period * 1.436 eg 17.23mo or 525 days
#    default=18 ; eg. 18 monthly periods : 50% attenuation at 15mo, close to zero pass at 12 mo. for 3 lobe.
#    lobes=3 ; use 525 for daily ; lobes=5 ; use 461 for 365 data
# console output comments give pass-band and stop band info. this is not sent to output file.

#    More lobes gives sharper transition band but more overshoot distortion / ringing.
#    default=3; 3 and 4 are good. ### >=5 may produce significant ringing artefacts ###

function sinc(x) {
  if (x==0) return 1; else return (sin(x)/x);

BEGIN{  OFMT = "%8.6f"

  if ( ARGC >3 ) {lobes=ARGV[3];ARGV[3]=""} else {lobes=3};
  if ( ARGC >2 ) {period=ARGV[2];ARGV[2]=""} else {period=18};
  if ( ARGC <2) { print "### usage: ./lanczos_LP.awk filename [50%_attn_period=18] [lobes=3]"; err=1; exit 1; }
  else {len=ARGV[2]; ARGC=2;}

  pi= 3.14159265359811668006 ;
  kwid= 2*halfk;         # kernel = -halfk....+halfk = kwid+1 pts ie odd number
  w_taper=twopi/kwid;    # first lanczos lobe -pi to +pi to fill window
  w=twopi/period;        # sinc passes zero pi, 2*pi, 3*pi ....

# calculate normalised kernel coeffs
  for (tot_wt=j=0;j<=halfk;j++) {tot_wt+=gwt[-j]=gwt[j]=sinc(j*w)*sinc(j*w_taper)};
  for (j=-halfk;j<=halfk;j++) {gwt[j]/=tot_wt};

# strip off last .xxx part of file name
# improve this  (doesn't work with paths like ../filename)

  if (ARGV[1]=="-") {
  else {
    for (i=1;i<length(fn);i++){   # if multiple dots in name, build back up except last part
      else basename=basename"."fn[i];
    print "# ",out_file >out_file;
    print "# ",out_file

  ln=-1;  # init linenumber counter

($0 !~ /^#/)&&($0 != "") {    # strip out comments, headers and blank lines

  if (ln>2*halfk)
    for (j=-2*halfk;j<=0;j++) {sum+=ydata[ln+j]*gwt[j+halfk]}
    if (out_file=="-")
      print xdata[ln-halfk],sum
    else {
      print NR,xdata[ln-halfk],sum
      print xdata[ln-halfk],sum >> out_file;
  } #endif ln


    if (err) exit err;
    if (FNR<=kwid) {print " ### insufficient data to fill "kwid+1" point filter buffer"; exit 1}
    print "# "lobes"-lobe lanczos kernel width = "2*halfk+1",done"
    if (lobes==5) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.258,period/.752)
    if (lobes==4) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.325,period/.698)
    if (lobes==3) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.436,period/.588)

    print "# output file = "out_file
# uncomment following to dump kernel (impluse response):
# sum=0; for (j=-halfk;j<=halfk;j++) {sum+= gwt[j]; printf ("%.4f\n", gwt[j])}
# print "# abs sum of coeffs = ",sum;


simple running-mean script ( as base for triple RM )

The following script will calculate a simple running mean. On its own this is a very poor filter, see accompanying article :


This script is intended for use by the triple-running-mean script available here:


Select code with mouse to copy elsewhere.

#!/bin/awk -f

# calculate running mean of width 'window'
# line number taken as x if file if file is a single column of data
# [usage] ./runmean.awk filename <window>
# default window=12
# strips empty and comment lines

# operates on files with two columns of numeric data
#  eg date as integer year or floating point decimal
# check whether data is continuous and evenly spaced !!

BEGIN{  OFMT = "%10.8f"
# ARGV[1]=fn; argv[0] is script name, hence ARGC>=1

  if (ARGV[1]=="") { print " ### error, no input file given." ; err=1; exit 1}

  if  ( ARGC >2 ) {window=ARGV[2];ARGV[2]=""} else {window=12};
  if (int(window)!=window) {window=int(window+0.5); window_adj=1; }

  print "# filtering "ARGV[1]" with running mean  window = ",window

  if (ARGV[1]=="-") {
   else {

    print "# ",out_file >out_file;
}  # end BEGIN

{ # body
if ( ($0 !~ /^#/) && ($0 !~ /^[" ",\t]*$/)  )  {  # strip empty and comment lines

  if ($2==""){  # for single column data files


  if (ln>=window)
    mid_x=(xdata[ln]+xdata[ln-(window-1)])/2;  # -1 to centre correctly
    sum-= ydata[ln-window+1];  # remove oldest
    if (out_file=="-")
      print mid_x,ymean;
    else {
      print NR,mid_x,ymean
      print mid_x,ymean >> out_file;
  } # endif
 } # non empty
} # end body

  if (err==1) { exit 1}
  if (window_adj==1) print "# warning : window rounded to nearest integer, window  = ",window
  print "# running mean window = "window",done"
  print "# output file = "out_file