A representative temperature record for New Zealand based on station data from 1853 onward is used in conjunction with four coupled climate models to investigate the causes of recent warming over this small midlatitude country. The observed variability over interannual and decadal time scales is simulated well by the models. The variability of simulated 50-yr trends is consistent with the very short observational record. For a simple detection analysis it is not possible to separate the observed 30- and 50-yr temperature trends from the distribution created by internal variability in the model control simulations. A pressure index that is representative of meridional flow (M1) is used to show that the models fail to simulate an observed trend to more southerly flows in the region. The strong relationship between interannual temperature variability and the M1 index in both the observations and the models is used to remove the influence of this circulation variability from the temperature records. Recent 50-yr trends in the residual temperature record cannot be explained by natural climate variations, but they are consistent with the combined climate response to anthropogenic greenhouse gas emissions, ozone depletion, and sulfate aerosols, demonstrating a significant human influence on New Zealand warming. This result highlights the effect of circulation variability on regional detection and attribution analyses. Such variability can either mask or accelerate human-induced warming in observed trends, underscoring the importance of determining the underlying forced trend, and the need to adequately capture regional circulation effects in climate models.
The anthropogenic component of surface temperature change has now been separated from changes resulting from internal climate variability and natural forcings for both the global mean temperature (e.g., Stone et al. 2007) and smaller continental regions of the earth (e.g., Stott 2003; Zwiers and Zhang 2003; Karoly and Braganza 2005b; Karoly et al. 2003). However, detection and attribution at the model gridbox scale is still a challenging exercise. This is principally due to the higher variability in regional temperatures because of local effects that, in turn, may not be well resolved by climate models (Hergerl et al. 2007). Despite this, Karoly and Stott (2006) established that for one model [the third climate configuration of the Met Office Unified Model (HadCM3)], the trend in central England temperatures (CET) since 1950 was significantly larger (at the 5% level) than trends resulting from natural internal (unforced) climate variability, but was consistent with the trends simulated by the model with anthropogenic forcings included. The CET extends back to before 1700, providing a long record of observed temperatures with small industrial influences, thus allowing for an evaluation of the model’s simulation of natural internal variability.
This paper attempts a similar analysis for the country of New Zealand, but also investigates the important role that circulation variability can play in such an analysis. With two main islands and a land area of about 270 000 km2, New Zealand is usually represented in current climate models by just a few grid boxes. It is well understood that, for New Zealand, the strongest circulation effect on temperature is whether the mean circulation is more northerly or southerly (Trenberth 1976). In this paper we first address whether climate models can simulate this temperature variability well, and then we consider whether the trends in New Zealand temperature can be successfully attributed to anthropogenic forcings if this dominant circulation effect on temperature is removed from both the observations and the climate models.
There are some studies that have taken similar approaches to detection and attribution analysis. On a global scale, Wu and Karoly (2007) considered 100-, 50-, and 30-yr trends up to 2002 calculated from the Hadley Centre Climate Research Unit global temperature dataset, version 3 (HadCRUT3v) at 5° latitude–longitude resolution. They found that in many areas of the earth the observed warming trend was larger than can be explained by internal variability at the 5° grid scale. They also considered the case where temperature variability resulting from the large-scale modes of the atmospheric circulation is removed first (e.g., the Arctic Oscillation). This resulted in a small reduction in the number of grid boxes in which trends were detectable, but this number was still significantly greater than that which might be expected from internal variability alone. The New Zealand landmass was one of the areas where trends were not separable from internal variability in either case.
For Australia, Karoly and Braganza (2005a) used the strong interannual relationship between temperature and rainfall in a number of regions to reduce the variability in annual mean temperature before successfully performing a detection and attribution analysis. Similarly, Douville (2006) showed that in the Sudan and Sahel regions the local temperature trends have been modulated by precipitation variability. In New Zealand any relationship between temperature and precipitation is much weaker. In fact, using the precipitation datasets for New Zealand, described in Tait et al. (2006), we were not able to find a significant relationship between annual mean temperature and precipitation averaged over the New Zealand landmass. More generally, Trenberth and Shea (2005) showed that the negative correlation between rainfall and temperature holds generally for continental areas at low latitudes, because dry conditions favor more sunshine and less evaporative cooling, while wet summers are cool. Over the ocean and at high latitudes the correlations were weaker and seasonally dependant.
Here we consider a large-scale circulation-based index that, unlike precipitation, is strongly correlated with New Zealand annual temperature variability, and we investigate the effect of removing this variability from a detection and attribution analysis. We begin by outlining and analyzing the best long-term record of New Zealand temperatures available and the four different climate models used for the analysis.
a. The New Zealand temperature series
The New Zealand temperature series (NZTS) was first outlined in Salinger (1980a), but has been subsequently updated and revised. The following seven locations were used to create the series: Auckland, Masterton, Wellington, Hokitika, Nelson, Christchurch, and Dunedin. The series consists of uninterrupted monthly mean temperatures from 1853 to 2006, though stations are included in the series only as they became available. All seven stations are included from 1907 onward, increasing the representativeness of the record from this point. The stations used in the record were carefully chosen not only for their record lengths, but also in order to be representative of the temperature variations observed across New Zealand resulting from the complex topography and varying exposure to atmospheric and oceanic circulations (Salinger 1980b). For each station the data have been carefully considered with regards to homogeneity, quality of the observations, environmental changes in the surroundings, and the different types of instrumentation used (see Rhoades and Salinger 1993 for details). In this study, annual means are calculated from the monthly data, but rather than for calendar years we calculate the annual mean from December of the preceding year through to November of the assigned year. Hence, the first year in the annual mean NZTS used here is 1854 rather than 1853. This was necessary for consistency because some of the model data were only available as seasonal means, with summer defined as consisting of the months of December–February, and so on.
b. Climate models
We have included simulations from four climate models, all of which were included in the fourth assessment report of the Intergovernmental Panel on Climate Change. These include the two Met Office models—HadCM3 (Johns et al. 2003; Tett et al. 2002) and Hadley Centre Global Environmental Model version 1 (HadGEM1; Johns et al. 2006; Stott et al. 2006)—and two National Center for Atmospheric Research (NCAR) models—the Parallel Climate Model (PCM; see Meehl et al. 2004; Washington et al. 2000) and Community Climate System Model, version 3 (CCSM3; see Meehl et al. 2006; Collins et al. 2006). Full details of the models and the simulations are provided in the indicated references. The choice of these models was determined by the ready availability of simulations of three different types. For the first, long control simulations are required in which external forcings are kept constant, and, for the case of greenhouse gases, kept at preindustrial levels (PICNTRL). In the second set of runs, the most important climate forcings of the twentieth century are incorporated, including observed and estimated changes in atmospheric greenhouse gases, ozone and sulfate aerosols, and natural external forcings, such as changes in solar irradiance and volcanic aerosol amounts (ALL). In the third set of simulations, only the natural external forcings are included, and these represent the simulated climate of the twentieth century had there been no anthropogenic influences (NAT). The models we analyze have a range of climate sensitivities and different assumptions about forcings. For example, the PCM includes the direct but not the indirect effect of sulfates, whereas HadGEM1 includes greenhouse gases, the direct and indirect effects of sulfate aerosols, and additional forcings, such as industrial black carbon, biomass-burning aerosols, and land use changes. The models also included different specifications of solar and volcanic forcings.
For each model a number of different initial condition ensembles were produced for each of the ALL and NAT simulation types. Table 1 details the number of ensembles for each model and the differing simulation lengths available. In total, there were 16 ALL simulations and 16 NAT simulations, while only the HadGEM1 simulations extend beyond 1999. The periods chosen for the PICNTRL runs excluded times at the beginning of the simulations where the models initially exhibited significant drifts in average New Zealand temperatures. These drifts might reasonably be assumed to be due to model spinup. To be comparable with the NZTS, which is a land-based temperature record, all New Zealand grid boxes in the model with land fractions of less than 50% were excluded from the analysis. The number of grid boxes remaining for each model is also listed in Table 1.
a. The effect of circulation anomalies on New Zealand interannual temperatures
Since 1900, New Zealand annual mean temperature has risen by an amount that is very similar to the global mean (see Fig. 1). For New Zealand most of the warming has occurred after 1946. However, in comparison to the global mean, the interannual and decadal fluctuations in temperature in the New Zealand region are significantly larger. As such, the possibility of attributing the temperature change, as has been done for the global mean, is much harder. In Fig. 1 it is also apparent that a number of abnormally cold years have been preceded by large equatorial volcanic eruptions, as was recognized by Salinger (1998). Interpretation of these events is complicated by the fact that strong El Niño events followed the eruptions of Mount Agung, El Chichón, and Mount Pinatubo, and these in turn have significant impacts on the New Zealand climate (Mullan 1995). Salinger (1998) showed that after removal of the ENSO signal there was still significant cooling following these eruptions, and this was due to the occurrence of more southerly flows over the New Zealand region.
More generally, the influence of southerly and northerly flow variability on New Zealand temperatures has long been recognized. Trenberth (1976) defined weather regimes P1 and P2 as the first two principal components of an empirical orthogonal function analysis of regional mean sea level pressure. Temperature changes over New Zealand were found to be principally explained by fluctuations in P2, which described the degree of southwest–northeast flow. P2 was found to be highly correlated (r = −0.84) with a station-based index called M1, defined as the pressure observed at the Chatham Islands subtracted from the pressure observed at Hobart. As such the index is a measure of the meridional flow between 147.5°E (Hobart) and 176°W (Chatham Islands) near 43°S, and is positive when the flow is toward the north. A number of similar indices were defined to describe zonal and meridional flow in the New Zealand region, but none showed such strong correlation to temperature as M1. As further evidence of the importance of the mean circulation on temperature, Salinger and Mullan (1999) found that in the period from 1930 to 1950 lower temperatures were associated with more south-to-southwest anomalous flow.
Figure 2 shows the NZTS, as for Fig. 1, but where the bold overlay is smoothed with a low-pass filter with eight terms. Overplotted is the M1 index anomaly from 1897 (the beginning of the Chathams record) through 2006, with the same smoothing as that applied to the NZTS. The figure suggests, qualitatively, that there is an association between the M1 and the NZTS: temperatures over New Zealand are warmer when the M1 anomaly is negative and cooler when it is positive. On physical grounds this is likely to be due to the advection of large-scale air masses originating from distinctly different regions of sea surface temperature; to the north of New Zealand there are warm subtropical waters, while to the south is the significantly colder Southern Ocean.
The M1 is characterized by a decreasing trend until around 1960 and an increasing trend subsequently. In conjunction with the warming temperature trend since 1960, this suggests that New Zealand has warmed despite a trend toward a greater frequency of southerly flows. To demonstrate the relationship statistically, we choose to perform a linear regression on the annual M1 and NZTS data. For this we use only the data from the period of 1907–46. This shortened time period simplifies the regression because we assume that during this time any contribution to the temperature record resulting from greenhouse gases is negligible; that is, the temperature record in this period is assumed to take the form
where T is the temperature at time t, βM1 the coefficient for M1 (the pressure index defined earlier), and ɛ is a residual series. We neglect data prior to 1907 because of the need for corrections to the M1 data prior to July 1906, and because this is the point at which all seven stations are included in the temperature record. We note that using the full record has no effect on either the slope of the regression or the correlation coefficient. The data and the regression are shown as a scatterplot in Fig. 3. For the period of 1907–46 (solid dots) there is a strong anticorrelation between temperature and M1 (r = −0.69) such that 47% of the variance in temperature is explained by the interannual variability in M1. The short duration of the record results in high errors for β of about 34% (given here as the 95% confidence intervals).
It is necessary in any linear regression analysis to establish that there is no autocorrelation in the residuals, because the ordinary least squares method inherently assumes that the error terms are independent. If this is not the case, the standard error estimates will be poor. Here we have tested for first-order autocorrelation using the Durbin–Watson test. For this test a value of 2 indicates zero autocorrelation, with a number less than 1 or greater than 3 indicating considerable autocorrelation. The test statistic returned in this regression is 1.57, which can be formally compared with a Durbin–Watson upper bound at the 5% significance level of 1.54. Thus, there is reasonable statistical evidence that for these data the residuals are not autocorrelated.
The last 40 yr (1967–2006) of M1 and temperature is also shown in Fig. 3 as open squares. Though the slopes are different for the two periods they agree within errors. However, there is a clear offset between the two linear best fits. This suggests a change in the background temperature of the New Zealand region over a large area, with no significant change in the north–south temperature gradient, and therefore no large change in the relationship between temperature and M1 variability (which would cause a change in the slope).
Figure 4 is a scatterplot of M1 versus temperature for the HadCM3 preindustrial control simulation. The slope is somewhat less and the correlation is reduced compared to the observations, suggesting a somewhat weaker relationship between temperature and the M1. This is the case for all the models, as listed in Table 2; however, as shown by the 95% confidence intervals derived from the standard error of the regression, it is not possible to say that the models are significantly different for the very short observational record. A reduced slope suggests that the models might be simulating weaker north–south gradients in the regional sea surface temperatures, which are likely to be responsible for the relationship between the M1 and New Zealand temperatures.
We now limit ourselves to considering the temperature record from 1947 onward in order to investigate 30- and 50-yr temperature trends under the presence of global warming. Here we assume that the observed temperature series is a function of M1 as derived from the earlier part of the record, as well as a linear temperature trend τ(t) (30 or 50 yr only), which approximates the effects of greenhouse gases and other anthropogenic forcings on New Zealand temperatures once the temperature variability associated with M1 has been removed; that is,
Because of the uncertainty in β this results in a range of possible temperature trends. The result of removing the M1 variability from the observational time series (1947–2006) is shown in Fig. 5. The variance of the record is reduced by 27%. In addition, because of the increasing M1 over the last 50 yr of the record it is apparent that a linear trend fitted to the period of 1950–99 is greater in the M1-removed residual temperature record. The 50-yr trend that is shown is the chosen example only because it is used later in the detection and attribution analysis, and the result is true of any 50-yr trend in this period. Most of the very cold and very warm years are also reduced in amplitude. The correlation between M1 and New Zealand temperature for each of the model control runs has also been used in an identical way to produce 16 transient model simulations with M1 variability removed.
Ideally, it would be possible to demonstrate that the trend in M1 over the last 50 yr is itself not anthropogenically induced, but rather results from internal variability. In this case we would simply be increasing the signal-to-noise ratio by removing internal variability from the analysis. In Fig. 6 M1 trends from the natural and anthropogenic transient ensemble members appear to cover about the same range, and hence give no suggestion that the observed M1 trend might be expected to be anthropogenically driven. However, it is also clear that the observed 50-yr trend in the M1 lies in the tail of the distribution of 50-yr trends from the model preindustrial control runs. Hence, it is likely that the models either do not capture the process well, which has led to the trend in the M1 in the New Zealand region, or that they significantly underestimate the natural variability in M1. Given that the models underestimate the interannual variability in the M1 (Table 2), the latter is a distinct possibility. This emphasizes the importance of being able to capture regional circulation effects in climate models before attempting detection and attribution analyses. Here we circumvent this by attempting to compare the models to the observations after the effect of the M1 has been removed through linear regression.
b. Internal and forced variability
For a detection and attribution analysis it is important to demonstrate that the models are able to capture the unforced and forced variability of the regional temperature observations. Here we use two separate methodologies to evaluate this. First, we calculate the standard deviations of interannual and decadal means for 100 yr of the NZTS from 1854 to 1953. Second, the standard deviation of 50-yr trends is calculated by using overlapping 50-yr windows that are 10 yr apart (i.e., from 100 yr there are six correlated 50-yr trends). For the models, all of the years available for each PICNTRL simulation are used. The standard deviations of interannual, decadal, and 50-yr trends are calculated in an identical way to the observations, but are sampled throughout the preindustrial controls by using 100-yr segments that advance by 50 yr. The variability of these individual estimates taken from the control run is used to estimate the 90% confidence intervals. For the 50-yr trends this method results in estimates of the standard deviations that have low biases for both the models and the observations, because there are only a limited number of independent estimates of 50-yr trends in a 100-yr sample. However, if the observational standard deviation lies within this confidence interval of the model deviations, then we have some evidence that the model can be said to have variability over these selected time scales, which is consistent with the observational record. In Fig. 7 it can be seen that for the NCAR models the observed interannual variability does not lie within the modeled 90% confidence interval. For longer time scales all of the models are consistent with the observations. For the 50-yr trends the very short observational record means that the results should be treated with caution. This comparison only assesses consistency, and there is a suggestion in the results that the models may have a somewhat lower average variability than that of the observations. The failure of the models to simulate sufficient M1 trends in the observations may be a possible cause of this. Given the very short observational record it is not possible to determine definitively whether the models are able to capture the real-world variability of 50-yr temperature trends.
While there is likely to be some climate response from increasing greenhouse gases toward the end of this period (1953), we believe that this represents a harder test of the model’s variability than if a longer preindustrial observed record were available; that is, the observational variability will be higher because of the inclusion of a component that may be anthropogenically influenced. This may partially explain the discrepancy with the interannual variability. Even with a longer observational record, such as in Karoly and Stott (2006), the observations contain natural forcings while the model preindustrial control does not, so any comparison is imperfect.
For the case in which M1 variability has been removed the observed temperature record is considerably shorter (1897–2006). In this case we have calculated interannual and decadal standard deviations from the 57 yr that are available from 1897 to 1953 in the same manner as for Fig. 7. However, this period is too short to make an estimate of the variability in 50-yr trends, and so we use instead 100 yr from 1897 to 1996 and compare this to the variability of the 50-yr trends from the same period of the twentieth-century transient simulations with all forcings included, but with M1 variability removed. In this case the confidence interval is simply given by the maximum and minimum of the ensemble members for each model. Because of the uncertainty in M1 removal the observed variability now has a 95% confidence interval. Comparison over this time period is a test of the model’s ability to simulate the correct variability of 50-yr trends with both internal and forced variability. This is similar to the approach taken in Hergerl et al. (2007), where power spectra are calculated from transient simulations as a way to assess simulated variability. In Fig. 8 the variability of observed interannual, decadal, and 50-yr trends all with M1 removed is consistent with the variability in all of the model preindustrial simulations with M1 removed. Taken together, Figs. 8 and 9 give strong confidence in all four models’ ability to correctly simulate the interannual and decadal variability in New Zealand temperatures. In addition, the short observational record cannot be used to reject the assumption that the models are able to capture 50-yr trend variability. This, in turn, is encouraging for the detection and attribution analysis that follows for 50- and 30-yr trends.
c. Detection and attribution
To assess the reasons for the recent trends in the NZTS, Fig. 9 compares the observed linear trends over 50- and 30-yr periods, ending in 1999, with the probability distribution from the combined model control runs. This was possible because all of the models exhibited similar variability in their control runs (Figs. 7 and 8). The end date for the trends of 1999 is chosen simply because this is the year for which three of the models’ simulations finish. The same plots are repeated for the case where the M1 variability is removed. The 50-yr trend is only detectable at the 5% level for the case where the M1 variability has been removed. Formally, the detection is done under the assumption that the preindustrial control simulations can be substituted for reality. The detection result is hence dependent on the assumption that the models are able to represent the variability of the real world both with, and without, M1 removed. In section 3b it was demonstrated that this hypothesis cannot be rejected on the basis of the limited observations, but the results must be considered with this caveat in mind. Specifically, detection is a one-sided significance test under the assumption that the PICTRL distribution is Gaussian. There is no attempt to account for the uncertainty in the M1 removal for this test because this is accounted for by the assumption that the models represent reality after M1 is removed. The absolute numbers are given in Table 3. In both cases the 30-yr trend cannot be separated from internal variability. Figure 9 also includes the linear trends calculated in the same way from the model transient simulations for NAT and ALL (16 of each). The 50-yr linear temperature trend with M1 variability removed can be detected and attributed to anthropogenic forcings, because the trend is consistent with the all-forcing transient simulations but is not consistent with the simulations with natural forcings only. Again, the specific ranges for the simulations are listed in Table 3. Because ranges are used to make this attribution statement the results must include the caveat that using a larger sample of models would result in a larger range for both the natural and anthropogenic simulations.
Figure 10 follows the same methodology as that of Fig. 9, but for trends finishing in 2006. In this case only HadGEM1 has transient simulations available. This gives a much-reduced sample of models, but a similar result is achieved. The 50-yr trend with M1 variability removed can be detected (5% level) and attributed to anthropogenic forcings. The warming in the last 7 yr up to and including 2006 has also been enough to now detect and attribute the 30-yr trend at the 5% confidence level. While the limited number of ensemble members must be acknowledged, Fig. 10 highlights the robustness of the 50-yr trend over a slightly different time span, and the importance of the increasing warming in recent years. It should also be noted that the detection of 50-yr trends at the 5% level is robust to whichever control distribution is used.
New Zealand surface air temperatures and the M1 circulation index have been shown to be highly correlated. The trend of increasing M1 (southerlies) since 1960 remains unexplained, because the climate models considered do not simulate such a trend and give little guidance on whether the trend might be due to internal variability or external forcing. With M1 variability included, all four climate models somewhat underestimate the temperature variability at the 50-yr time scale when their preindustrial controls were compared to observations prior to 1954 (though all were still within the 90% confidence intervals). They also somewhat underestimate M1 variability, and consequently the model and observations are in better agreement if this variability is removed from both.
The presence of a positive trend in M1 since 1960 suggests that temperatures in New Zealand would have otherwise been warmer than observed, and may in fact increase rapidly if the trend were to stop or decrease. Both 50- and 30-yr temperature trends in New Zealand can only be reliably detected and identified as being consistent with anthropogenic forcings, but not consistent with only natural forcings, if M1 variability is first removed. Thirty-year trends ending in 1999 still cannot be separated from internal variability even after M1 variability has been removed. Thirty-year trends are more sensitive to endpoints, and were clearly increased when the M1 variability was removed.
There are a number of very strong caveats to these findings. The detection step relies on the assumption that the models adequately simulate variability on the regional scale. It was demonstrated that this was a reasonable assumption for interannual and decadal means both with and without M1 removed. The observed record was also found to be not inconsistent with the simulated 50-yr trend variability. However, it was also demonstrated that the models do not include a trend in the M1 as seen in observations, and this is evidence that on the regional scale this assumption may be poor. The analysis is improved by a simple approach to removing the M1 variability from both the observations and the models, but this approach cannot rule out the possibility of compensating errors. When a model failure is identified the best approach is to establish the cause of the possible problem and work toward improving the model simulations. The attribution step makes use of the ranges in trends from a limited number of twentieth-century simulations from four different models up to 1999, and only one up to 2006. Increasing the number of models/simulations would likely result in an increase in those ranges.
The analysis did not remove the effect of the southern annular mode (SAM) on New Zealand temperatures. Gillett et al. (2006) showed that some station temperatures in New Zealand exhibit a correlation with the phase of the SAM that is robust at the 90% significance level. Wu and Karoly (2007) presented results that indicated a reduction in the warming in the New Zealand region when the effects of the SAM, AO, and Pacific–North America (PNA) pattern are removed from 30- and 50-yr temperature trends derived from gridded temperature records. However, for annual mean temperatures for the combined station data used in this work there is only a small positive correlation (r = 0.22) between the index of the SAM and temperature. This effect is significantly smaller than the M1 correlation but could have been included in the analysis. However, its interpretation is complicated by the fact that the trend in the SAM is likely to be associated with anthropogenically driven forcings caused by ozone depletion and greenhouse gas increases (Arblaster and Meehl 2006).
This study supports the hypothesis that a regional detection and attribution analysis may be more successful when local circulation effects, which are poorly captured by models, have been removed. The simple, consistency type of attribution analysis done here could be improved by undertaking a full optimal detection analysis incorporating an M1 pattern, as Gillett et al. (2000) employed for analyzing the effect of the North Atlantic Oscillation. Identification of the cause of the M1 trend will be a priority for further work, particularly with regard to understanding any possible effects on future climate change over New Zealand.
The authors thank Brett Mullan and James Renwick for discussion on the ideas developed in this paper, and Jonathan Gregory for assistance in acquiring some of the climate model data. We also acknowledge the significant improvements contributed to this paper by the anonymous reviewers. This work was partially funded by the New Zealand Foundation for Research Science and Technology under Contract C01X0701 and by the Defra and MoD Integrated Climate Programme GA01101, CBC/2B/0417_Annex C5.
Corresponding author address: Sam Dean, Private Bag 14901, National Institute of Water and Atmospheric Research, Kilbirnie, Wellington, New Zealand. Email: email@example.com