Updated: added residuals plot

I’m really impressed by this. I now have a snow depth prediction model that “explains” 50% of the variance in the Spencers Creek season peak snow depth record (midway between Perisher Valley and Thredbo, from Snowy Hydro). When you consider the vagaries of weather, snowfall, compaction, melt and weekly measurement, that really is a surprising achievement. I sure didn’t expect to get anywhere near it when I started out nearly a decade ago. Here’s how the new model goes compared with last year’s effort (mark III):

Like before, the model is based on multiple linear regression with well known climatic variation modes and influences, as listed below. As promised, I’ve added the relationship with southern hemisphere volcanic aerosols, but eschew the dubious (and rather weak) correlation with sunspots. You can see the improvement in the big volcano years: Agung in 1964 and Pinatubo in 1991 and 1992.

I’ve also extended the regression “training” period from the previous start in 1958 back to 1955, using IOD estimates for those years from Cai et al (2009). I still leave out 1954, because I’ve become convinced that the sparse depth measurements in that first year of our snow depth record — nearer monthly than weekly — may have badly missed the peak. (By far the worst model miss-fit in 61 years of data is in 1954: by 1.65 m, or nearly four standard deviations.)

The individual correlations look like this:

Despite the scatter, each of those is statistically significant at least at the 5% level (single tailed, because in every case the trend direction was clear a priori). Most are highly significant (p < 1%).
The prediction model equation is:
**Spencers Creek peak depth (cm) = 899 – 0.311 x year – 13.8 x AAO + 2.09 x SOI – 8.80 x IOD – 5.36 x PDO – 143 x SST + 534 x aerosol**

R-squared is 0.50 and the standard error is 47 cm. The parameters are as follows:

**Calendar year**, to model the downtrend. The calendar year term does not encompass all of the observed snow depth downtrend because there are substantial trends in other parameters, for example IOD (see the reference below) and AAO (see AR5 WG1, chapter 14.5.2) — although note that much of the relatively strong recent trend in the latter is usually attributed to ozone depletion (but mostly in summer). The even stronger uptrend in SST does not affect the model because I use a detrended SST series (see below).

**Antarctic oscillation**(AAO), also called “southern annular mode” or SAM, is a measure of how tightly the circumpolar winds (“polar vortex” in one usage) blow around the pole. A loose pattern (negative AAO) leads to more polar storms reaching southern Australia, and more snow. The winter average AAO is used — the average of June, July and August.

**Southern oscillation index**(SOI) is the difference between Tahiti and Darwin surface atmospheric pressure expressed as monthly standard deviations times ten. SOI is an indicator of the El Niño Southern Oscillation (ENSO), an east-west quasicycle in equatorial Pacific Ocean surface temperature and wind patterns which correlates with precipitation across much of Australia, including with alpine snow. A positive SOI is associated with more (and wetter) Australian snow. The winter average is used.

**Indian ocean dipole**(IOD) is an ENSO-like variation in the smaller Indian Ocean, which correlates with precipitation across southern Australia, including with alpine snow. Negative IOD is associated with more snow (the sign is consistent with SOI, but we’re on the opposite side of the Indian Ocean). The winter average is used.

**Pacific decadal oscillation**(PDO) is a long-cycle, mostly north-south variation in the western Pacific Ocean. Negative long-average PDO is weakly correlated with more snow. The 2-year average to August is used.

**Sea surface temperature**(SST) is that in the Great Australian Bight and northwest Tasman Sea, averaged over the box: latitude 30-37°S, longitude 115-160°E, and expressed as degrees Celsius anomalies from the 1951-1980 mean,*detrended*about 2015. Cool SSTs correlate strongly with more snow, but unfortunately local SSTs are rising rapidly with global warming. The winter average is used.(I detrend the SSTs used in the model to make the equation parameters appear more sensible, and to make them more comparable with previous versions. Because of the model linearity, linear detrending doesn’t alter the regression outcome — it just shifts the SST trend effect across to the “calendar year” term.)

**Aerosol optical depth**is that for the southern hemisphere at 550 nm. Large optical depths correlate with big snow seasons, while smaller depths appear to have little effect. The winter average is used.

### Residuals

Model residuals are the bits the model doesn’t get — the differences between the real data and the predictions. Those look like this:

The one standard deviation error is 44 cm, a little less than the reported standard error of the regression (the definitions are slightly different). It’s interesting that there appears to be no skew at all in the residuals, despite a modest level of skewness in the raw data. It’s possible that much of the raw data skewness derives from the occasional big volcanic years (which are now represented in the model). The oddly rectangular shape is also interesting (nearly all the errors fall withing ±70 cm — the bar labels are the midpoints), but that’s probably just a chance artifact.

To make my 2015 prediction, all that remains is to estimate the parameters and plug them in. As regular readers will know, that is stuff for next month.

### Reference

Cai, W., Cowan, T., & Sullivan, A. (2009). Recent unprecedented skewness towards positive Indian Ocean Dipole occurrences and its impact on Australian rainfall. *Geophysical Research Letters*, 36(11).