Section 1: Exploratory Data Analysis, Aggregation, and Decomposition of Sea Ice Volume

Exploratory Data Analysis

#Read in the sea ice volume data from UW's PIOMAS project
seaice<-readRDS("piomas.rds")

#Checking the parameters of the sea ice data to see start, end, and frequency. This data set was taken at a monthly frequency from 1979 to 2016.
tsp(seaice)
## [1] 1979.000 2016.917   12.000
ts.plot(seaice, ylab=expression(Seaice~10^3~km^3))
title("Sea Ice by Volume Measured by PIOMAS")

Fig 1: Arctic sea ice volume in thousands of kilometers cubed from PIOMAS project data every month from Jan 1979 to Dec 2016

summary(seaice)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.787  14.893  20.163  19.894  25.524  32.951

The average sea ice volume is 19,890 cubic km , with a median of 20,160 km which are fairly close, with the minimum sea ice volume observed at 3,787 cubic kilometers

Lets make a Normal Q-Q plot and see if the data is normally distributed

qqnorm(seaice,main="Normal Q-Q of Seaice Volume from PIOMAS")

Figure 2: Normal Q-Q plot of Arctic sea ice volume to look for a normal ditribution The sea ice data is fairly normally distributed, sample and theoretical quantiles match, forming a straight line and there is little skewness.

Next lets run a linear model to see what kind of linear trend we are dealing with. I am expecting the sea ice volume to show a negative slope overall.

seaicetime<-time(seaice)
seaice.lm <- lm(seaice~ seaicetime)
summary(seaice.lm)
## 
## Call:
## lm(formula = seaice ~ seaicetime)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1812  -5.9318   0.2425   5.6234   9.7438 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 644.74933   51.51744   12.52   <2e-16 ***
## seaicetime   -0.31275    0.02578  -12.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.04 on 454 degrees of freedom
## Multiple R-squared:  0.2447, Adjusted R-squared:  0.2431 
## F-statistic: 147.1 on 1 and 454 DF,  p-value: < 2.2e-16
ts.plot(seaice, ylab=expression(Seaicevolume~10^3~km^3), xlab="time")
title("Sea Ice by Volume Measured by PIOMAS")
abline(seaice.lm, col="red",lwd=2, lty="dashed")

Figure 3: Sea ice volume measured monthly from Jan 1979 to December 2016 with a linear model trendline showing a overall decrease in sea ice volume with time. The slope showed a decrease of 312.7 cubic kilometers per year since 1979.

Next I made a correlogram of the autocorrelation

acf(seaice, main="Autocorrelation in monthly Sea Ice Volume 1979-2016")

Figure 4: The correlogram of sea ice volume measured monthly from Jan 1979-Dec 2016, showing positive and negative autocorrelation which indicates a seasonal cycle.

Next lets take a look at the residuals

seaice.residuals<-residuals(seaice.lm)
seaice.residuals<-ts(seaice.residuals)
plot(seaice.residuals, main="Residuals of Sea Ice Volume by Month 9179-2016")

Figure 5: Plot of residuals of sea ice by volume in thousands of cubic kiolmeters from 1979-2016.

It appears that the Arctic sea ice volume is decreasing by 312.7 cubic kilometers per year over the last 38 years. The data are relatively normally distributed and show little skewness but there is a seasonal cycle of sea ice volume minimum and maximum every year. Since there is an overall negative trend in the data, it is not stationary so detrending or differencing will be needed later to interpret it. The correlogram (Fig 4) shows a significant positive autocorrelation followed by significant negative autocorrelation with great persistence or an “echo”" in at lease the first two years.In other words, the seaice data are positively autocorrelated for a few months and then negatively correlated for a few months and (hot or cool) tend to follow like months. The residuals show very little shape but toward the last 100 months the range of residuals is slightly larger, indicating more variation in the last 8 years or so (Fig 5).

Aggregation

Now it is time to aggregate our montly data into three different sets of annual data

yrs <- 1979:2016
# This is 12-month filter 
seaice.12.filt <- filter(x = seaice,filter = rep(x=1/12,times=12),sides=2)
head(cbind(seaice,seaice.12.filt))
##      seaice seaice.12.filt
## [1,] 27.704             NA
## [2,] 30.171             NA
## [3,] 32.045             NA
## [4,] 32.951             NA
## [5,] 32.295             NA
## [6,] 29.787       25.42608
#now its time to extract the average from the proper month, since we used sides=2 that would be halfway between OCt and Sep, or month March

iceyear.index <-seq(from=3,to=length(seaice.12.filt),by=12)
iceyear <- seaice.12.filt[iceyear.index]
#Length of iceyear object should be 38
length(iceyear)
## [1] 38
#Success!
##convert this object into a ts 
iceyear<-ts(iceyear)
iceyeartime<-time(iceyear)
lm.iceyear<-lm(iceyear~iceyeartime)
#Lets pull out the minimum of each year, September values
#I found a nifty package to subset either seasons, quarters, or months and retain the ts designation

library("forecast")
annual.minimum<-subset.ts(seaice,month=9)
annual.minimum.time<-time(annual.minimum)
lm.annual.minimum<-lm(annual.minimum~annual.minimum.time)
#Lets get the annual maximum from March using the same nifty package
annual.maximum<-subset.ts(seaice, month=3)
annual.maximum.time<-time(annual.maximum)
lm.annual.maximum<-lm(annual.maximum~annual.maximum.time)
#Lets make a cohesive plot
plot(iceyear, ylim=range(2,38),ylab="",col="green",lwd=2.5)
abline(lm.iceyear,lty="dashed",lwd="2",col="green")
par(new=TRUE)
plot(annual.minimum, ylim=range(2,38), axes = FALSE, ylab="",col="red",lwd=2.5)
abline(lm.annual.minimum,lty="dashed",lwd="2",col="red")
par(new=TRUE)
plot(annual.maximum,ylim=range(2,38),axes=FALSE,main="Seaice volume 1979-2016",ylab=expression(Seaicevolume~10^3~km^3),col="blue",lwd=2.5)
abline(lm.annual.maximum,lty="dashed",lwd="2",col="blue")
legend(38,35,legend=c("iceyear","annual minimum", "annual maximum"), col = c("green","red","blue"))

Figure 6:Ice volume over an ice year from Oct to Sep, annual minimum, and annual maximum from Jan 1979-Dec 2016. Data courtesy of UW’s PIOMAS project on Arctic sea ice.

The annual ice year values ranged from 14,000-24,000 cubic kilometers which was a range in between the annual minimum ice volume and the annual maximum ice volume. The plot of ice year seems like a good average of the two and all three series show similar negative trends. There are peaks for minimums in 1983, 1996, and 2012 and a positive peak over the trendline in 2015.The Annual minimum dipped farther below the trendline after 2010, indicating some mechanism may be at play making the minimum even lower than “normal”. Overall each data set had a similar shaped and similar trend but was just different in magnitude on the y axis.

Decomposition

Decomposition following an additive time series model

seaice.d <- decompose(seaice) 
plot(seaice.d)

Figure 7: Decomposition of Arctic sea ice volume from 1979-2016 showing the trend, seasonal component, and remainder “noise”.

The m-hat trend shows an overall negative trend in Arctic sea ice volume, reaching its lowest point just after 2010 and 206. The seasonal component corresponds to the cycles observed in the original time series, with minima in September and maxima in March. The residual noise “z-hat” covers a much larger ranch starting in 2010, indicating an increase in variation after 2010.

Deseasonalizing Arctic ice volume

deseasonalized<-aggregate(seaice,FUN=mean)

plot(deseasonalized, yax=T,ylab="Arctic Sea Ice Volume in 10^3km^3", xlab="time", main="Deseasonlized Arctic Sea Ice Volume")

Fig 8: Deseasonalized Arctic Sea Ice Volume from 1979-2016. Data courtesy of UW’s PIOMAS project.

The deseasonalized plot shows an overall downward trend and again minimum peaks at 1983, 1996, and just after 2010 and 2016. It is easier to decifer the overall trend and anomalies without the seasonal cycles in there.

Stationarity

Stationarity

Zhat residuals show icnreased range post 2005, which means that the time series is not stationary because the residuals are not constant throughout the time series, they are starting to change in magnitude. According to Duarte et al. 2012 and Tilling et al. 2015, the lowest 5 minima for Arctic sea ice volume occurred from 2007-2011. In pre-2007 years the ice sheets did not melt completely over the summer, but as summer ice volume drops lower and lower, the amount of ice for the cool season to add on to also drops (Duarte et al. 2012). With multi-year ice bases getting reduced to seasonal ice sheets, the variation in sea ice volume increases from year to year (Duarte et al. 2012). However, this may not be a “tipping point” and could be reversible according to studies carried out by Tilling et al. In 2013. Tilling found that there was a 25% and 33% increase in Arctic sea ice volume in the Autumn of 2013 and 2014 (particularly cold years), which reduced the number of days where melting occurred in the following summer (Tilling et al. 2015). This may offer a glimmer of hope if the time series ended in 2015, however, by the minima in 2016 it appears the Arctic sea ice volume has not continued to increase. In fact, the variability in sea ice volume has increased, making forecasting more unreliable.

Part 2:Autocorrelation and ARMA models

Background Reading

LaMarche 1974

LaMarche found that the high frequency patterns in the ring width index at all the sites were similar, showing a period of less than 10 years. The low frequency pattern present in the upper forest boundary was on the order of 27 years. Tree ring indices can provide an excellent proxy for studying the climate from thousands of years ago. The Snake Range of Eastern Nevada demonstrates a change in precipitation and vegetation with altitude, with Bristlecone Pine confined to the subalpine region. Sites were chosen to get different altitudes and slopes, with cores taken over the course of two years, cross-dated, then the ring width was measured. Standard deviation, mean sensitivity, and correlation coefficients were calculated for samples from all sites and revealed that the ring width index at the lower and upper forest borders, which had the most climatic variation, showed the largest response. Digital filters allow us to highlight patterns at different frequencies so the sites can be easily compared to analyze long term and short term patterns. High pass filters remove low frequency, long term patterns like year to year growth and save the high frequency peaks. The low pass filters remove the higher frequency patterns, leaving you with long term trends. LaMarche used both low and high pass filters on the chronologies from all sites and compared them with linear product-moment correlation coefficients. Comparing the high pass filtered time series of the four sites (including lower forest border and upper forest border) found similar high frequencies patterns amongst all four. However, looking at the low pass filtered time series differences between low forest border and upper forest border sites became apparent. All the lower forest sites had similar long term variation which correlated in a slightly negative fashion to upper forest sites. This means there is something going on at a low frequency that makes the RWI respond differently based on altitude. Principal components analysis was then used to evaluate individual trees at each site. PCA of the high pass filtered time series supported the correlation analysis and found similarities in high frequency variation among individual trees and among sites. Individual trees from each site behaved similarly to each other, indicating site was a good predictor of RWI response depending on the frequency analyzed. Spectral analysis can help identify periods in the data rather than fitting an arbitrary filter. Often, you will see an “echo” of a time series persist at the lower frequency end. LaMarche ran a spectral analysis which revealed the lower forest border chronologies tended to reach maxima or minima about three years before the upper forest border chronologies. Next mean temperature and precipitation data from a nearby site were used to model the RWI response to climate. LaMarche surmised that the change in tree ring character at the turn of the 19th century could be due to increased water use efficiency correlated to higher atmospheric carbon dioxide at the dawn of the industrial age. A well-known pattern of variation in wind and pressure data known as the “biennial pulse” cycles every 2.2 years which could account for the high frequency variation, indicating the lower forest border responds to changes in precipitation.

Hughes and Funkhouser 2003

Hughes and Funkhouser wanted to continue part of LaMarche’s work and expand it to other regions to see if the same patterns held over a larger range. They found higher autocorrelation for trees in the ufb versus those in the lower forest border, which is similar to what LaMarche foudn. Each mountain pair showed similar 3-7 year series, however, this author found no negative correlation between ufb and lfb groups when looking at lower frequency, multidecadal time scales. The 3-7 year series was significantly correlated with precipitations patterns, leading the authors to conclude that the lfb tree rings growth is limited by reduced precipitation. This is similar to what LaMarche found, however these authors had access to higher resolution climate data which bolstered their claim. UFB trees across sites showed similar low frequency patterns. The authors agree with LaMarche that this may be indicative of a temperature signal or the rise in carbon dioxide combined with increased temperature since the late 19th century. In conclusion, they found a shared low frequency precipitation signal among lfb and ufb and ufb alone showed a low frequency pattern corresponding to temperature/ carbon dioxide.

Plotting the two time series

lfb<-readRDS("lfb.rds")
ufb<-readRDS("ufb.rds")
par(mfcol=c(2,1))
par(mar=c(2,2,2,2))
ts.plot(ufb, ylab="RWI Upper Forest Border",yax=T, main="RWI for Lower and Upper Forest Border of BCP", xlab="year")
ts.plot(lfb,yax=TRUE, xlab="year",ylab="RWI lower")

Fig 9: Time series of RWI for upper and lower forest border of Bristle Cone Pine 1000-2009 CE

Autocorrelation and partial autocorrelation

acf(lfb,main="Autocorrelation of Lower Forest Border BCP")

Fig 10: Autocorrelation of lower forest border BCP 1000-2009CE

pacf(lfb, main="Partial Autocorrelation of Lower Forest Border BCP")

Fig 11: Partial Autocorrelation of Lower Forest Border BCP 1000-2009CE

acf(ufb, main="Autocorrelation of Upper Forest Border BCP")

Fig 12: Autocorrelation of Upper Forest Border BCP 1000-2009CE

pacf(ufb,main="Partial Autocorrelation of Upper Forest Border BCP")

Fig 13: Partial Autocorrelation of Upper Forest Border BCP 1000-2009CE

The ACF of the Lower Forest Border shows a value of 1 for the lag 0 autocorrelation so we may compare other lags to the maximum phi value of 1 (Cowpertwait and Metcalf 2008). Using the 95% confidence interval we can see that there is only significant autocorrelation at lag 1 and mildly at lag 2 according to the 95% confidence interval, with a cosine shape of shallow amplitude. When combined with the PACF plot showing signigicance at lag 2 this series has the hallmarks of a second order autoregressive model. The PACF in figure 11 shows that even after accounting for autocorrelation due to a first order model, there is still slightly significant autocorrelation, which led me to the conclusion of a second order model. The ACF of the Upper forest border (fig 12) looked very different from that of the lower forest border. There was a high serial correlation that gradually decayed overtime, but the PACF (fig 13) did not show a fall to 0 after lag 1, leading me to rule out the possibility of a random walk. Ecologically speaking a random walk didn’t make sense to me either for tree rings. In fact, the PACF (fig 13) shows an almost geometric decay, leading me to believe that we have an AR (1) MA (1) process at play or an AR 5 model because the PACF drops to 0 after lag 5.

ar(lfb)
## 
## Call:
## ar(x = lfb)
## 
## Coefficients:
##      1       2  
## 0.1370  0.0719  
## 
## Order selected 2  sigma^2 estimated as  0.0598

The AR function for the lower forest border confirms my interpretation from the correlograms that we are dealing with a second order model for the LFB with coefficients 0.1370 and 0.0719. The UFB is a little more difficult.

x.arma<-arima(ufb,order=c(5,0,0)) #A 5th order AR model
AIC(x.arma)
## [1] -997.0653
x.arma2<-arima(ufb,order=c(5,0,1)) # An AR(5) MA(1) model
AIC(x.arma2)
## [1] -999.9682
#This didn't make for a better AIC and made the model more complex, lets try a lower order

x.arma3<-arima(ufb,order=c(1,0,1))
AIC(x.arma3)
## [1] -1002.134
#This did not decrease our AIC or fit our correlogram
#How about trying an MA model only? 
x.arma4<-arima(ufb,order=c(0,0,3))
AIC(x.arma4)
## [1] -728.4874
#Coefficients:
    # 1       2       3       4       5  
#0.3036  0.2096  0.1253  0.1656  0.0626

Based on relative AICs and fitting multiple AR and ARMA models, the most parsimonious model was an AR(2) process for the Lower Forest Border and an AR(5) for the Upper Forest Border. The correlograms for the UFB clearly show strong autocorrelation, stronger than that taking place on the lfb so tree ring growth at time step 1 will affect tree ring growth five time steps in the future. These results are consistent with what LaMarch and Hughes and Funkhouser discovered in that a higher autocorrelation exists at the ufb sites and there are similar high frequency pattterns in both data sets. This high frequency pattern is likely due to the “biennial pulse” (LaMarche 1974) altering precipitation patterns and limiting lfb growth. This would explain why the lfb does not show as high autocorrelation as the ufb, as precipitation changes every 2-3 years so the RWI at time step t would be a function of all the steps until t-3 only, in the absence of some external MA process which I did not uncover.

In contrast, the ufb had much higher autocorrelation whereby the RWI five time steps away would still be affected t1 likely due to the slower process of temperature and carbon dioxide change. Schulze et al. (1967) provided experimental evidence that the bristlecone pine photosynthesis and growth is halted at lower temperatures, noting that over summer a high daily photosynthesis rate must be achieved to offset a negative CO2 balance incurred over winter and continue to grow.This explains why the UFB trees had a higher serial autocorrelation than the lower forest border because the tree’s carbon dioxide and temperature experience over winter can affect the RWI multiple time steps in the future. LaMarche (1984) further showed that bristlecone pine had enhanced growth with increasing temperature until the 1960s, whereby a continual increase in RWI was attributed to carbon dioxide fertilization.

Section 3: Phenology Regression of Cherry Blossoms

cherry<-readRDS("cherry.rds")
plot(cherry, main="Phenology of Cherry Blossoming")

Fig 14: Time series plot of Day of Year, March Temperature, and March Precipitation for Cherry Blossoms in Kyoto, Japan for 500 years

acf(cherry)

Fig 15: Autocorrelation of the three time series, Day of Year, March Temp, and March Precipitation along with the cross correlations

pacf(cherry)

Figure 16: Partial autocorrelation of three time series, Day of Year, March Precipitation, and March Temperature, for Tokyo Japan along with cross correlations

Looking at the cross correlation of Day of Year and March Temperature is interesting in that you can see at negative lags, March temperature before time T is correlated with Day of year at time T. I am interpreting this as the March Temperature before the day of blossoming helps predict, or leads, date of blossoming. The Day of year has a significant autocorrelation at lag which is to be expected of biological systems and March Precipitation was similar. March temperature was highly autocorrelated for multiple lags, indicating the temperature measurements at time T were a large function of the temperature at T-1, moreso than day of year and precipitation.

pacf(cherry)

Figure 17: Partial autocorrelation of Day of year, March Precipitation and Temperature of Kyoto, Japan

The only partial autocorrelation with significant lags was Day of Year and March Temp partial cross correlation and I am not sure how to interpret this.

Next I will separate the multiple time series into three ts objects.

doy.ts<-ts(cherry[,1],start=1882,freq=1)
marcht.ts<-ts(cherry[,2],start=1882,freq=1)
marchp.ts<-ts(cherry[,3],start=1882,freq=1)
par(mfcol=c(3,1))
par(mar=c(2,3,3,2))
acf(doy.ts, main="ACF for Day of Year Cherry Blossoms Bloom")
acf(marcht.ts, main="ACF for March Temp Kyoto, Japan")
acf(marchp.ts,main="ACF for March Precipitation, Kyoto, Japan")

Figure 18: Autocorrelation for Day of Year, March Temperature, and March Precipitation in Kyoto, Japan

par(mfcol=c(3,1))
par(mar=c(2,3,3,2))
pacf(doy.ts, main="PACF for Day of Year Cherry Blossoms Bloom")
pacf(marcht.ts, main="PACF for March Temp Kyoto, Japan")
pacf(marchp.ts,main="PACF for March Precipitation, Kyoto, Japan")

Figure 19: Partial autocorrelation of Day of Year, March Temp, and March Precipitation in Kyoto, Japan

The March temp had a higher serial correlation than the other two variables, which showed little significant autocorrelation after lag 1.

First I will run an ordinary least squares multiple regression with day of year for blossoming as a function of March precipitation and March temperature

ols1<-lm(doy.ts~marcht.ts+marchp.ts)
summary(ols1)
## 
## Call:
## lm(formula = doy.ts ~ marcht.ts + marchp.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.2013 -2.4502 -0.1654  1.9944  9.5266 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 125.4085    20.3195   6.172 1.13e-08 ***
## marcht.ts    -9.3876     1.0018  -9.371 1.00e-15 ***
## marchp.ts     0.6071     0.3395   1.788   0.0765 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.796 on 111 degrees of freedom
## Multiple R-squared:  0.4419, Adjusted R-squared:  0.4318 
## F-statistic: 43.94 on 2 and 111 DF,  p-value: 8.767e-15

Next I will look at the ACF and PACF plots of the residuals for the ols1 model to see if there is autocorrelation in the residuals. If there is, then the residuals are not independent on one another than this assumption is violated.

resi<-residuals(ols1)
acf(resi)

pacf(resi)

auto.arima(resi)
## Series: resi 
## ARIMA(2,0,2) with zero mean     
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.0124  -0.9597  0.0679  0.8623
## s.e.   0.0427   0.0764  0.0871  0.1487
## 
## sigma^2 estimated as 12.53:  log likelihood=-304.31
## AIC=618.62   AICc=619.18   BIC=632.3

The ACF plot in fig ? shows a cycle of significant positive and negative autocorrelations which makes it appear there is not just an AR process but also an MA (external) process. The PACF shows significant autocorrelation after lag 1 and doesn’t fall to 0 so we have something more than an AR 1 process at work. Even after accounting for autocorrelation at the previous lag, we still have significant autocorrelation at multple lags. Next, I used auto.arima to fit an ARMA model to the ols1 multiple regression model. I found that the best fit was an AR2, MA 2 process.

Section 4: Ocean Patterns and Spectral Analysis, Smoothing

First I will read in the data and split it into separate time series

oceans<-readRDS("OceanPatterns.rds")
plot(oceans, main="ocean circulation patterns")

Figure 20: Ocean Patterns for ENSO, PDO, and AMO

tsp (oceans)
## [1] 1900.000 2016.917   12.000
enso.ts<-ts(oceans[,1],start=1990,freq=12)
pdo.ts<-ts(oceans[,2],start=1990,freq=12)
amo<-ts(oceans[,3],start=1990,freq=12)
##separating the ts from mts

Applying a Lowess smoother

enso.lo<-lowess(enso.ts,f=1/15)
plot(enso.ts, yax=T, main="ENSO")
lines(enso.lo,col="red",lwd=2)

Figure 21: El Nino Southern Oscillation with lowess smoother

pdo.lo<-lowess(pdo.ts,f=1/10)
plot(pdo.ts,yax=T,main="PDO")
lines(pdo.lo,col="green",lwd=2)

Figure 22: Pacific Decadal Oscillation with Lowess smoother

amo.lo<-lowess(amo,f=1/5)
plot(amo,yax=T, main="AMO")
lines(amo.lo,col="blue",lwd=2)

Figure 23: Atlantic Multidecadal Oscillation with Lowess Smoother

I decided on different smoothers based on the more low frequency patterns I observed. I noticed that the Enso had a noticeable pattern that was at a higher frequency than the other two. The PDO was next, followed by the lowest frequency pattern noticeable in the AMO. The lowess filter has the benefit of keeping the end points, unlike a moving average or a Hanning, which is useful in this case because these are relatively short time series. The f factor chosen is the fraction of points used to calculate a polynomial fit, so I chose smaller denominators for the processes I thought contained lower frequency patterns.

enso.spec <- spectrum(enso.ts, span=5, plot= FALSE)
plot(enso.spec,log="no", type="h", xlab="Frequency (cycles / yr)", ylab="Spectral density",main="ENSO Ocean Patterns",xlim=c(0,2.05))

#annual variation?

Figure 24: Periodogram of ENSO showing peaks in cycles per year

pdo.spec<-spectrum(pdo.ts,span=5,plot=F)
plot(pdo.spec,log="no",type="h", xlab="Frequency (cycles/yr)",ylab="Spectral density", main="PDO Ocean Patterns",xlim=c(0,0.6))

Figure 25: Periodogram of Pacific Decadal Oscillation

amo.spec<-spectrum(amo,span=5,plot=F)
plot(amo.spec,log="no",type="h", xlab="Frequency (cycles/yr)",ylab="Spectral density", main="AMO Ocean Patterns",xlim=c(0,0.6))

Figure 26: Periodogram of Atlantic Multidecadal Oscillation

The ENSO showed a dominant frequency of 1 which corresponds to a 1 year period. This implies an annual pattern for ENSO and although there are smaller peaks present at lower frequencies, these are likely a harmonic remnant of the dominant frequency at 1.The PDO periodogram showed dominant frequencies at 0.02 and 0.03 which correspond to periods of 30-50 years. There are higher frequencies at 0.1 and 0.2 which correspond to 3 and 5 year periods, though they are smaller than the frequency peak at 0.02. The AMO periodogram had a dominant frequency at 0.01 through 0.03 that were by large greater than any other frequencies. This corresponds to a period of 50-100 years which is a much longer period than the other two patterns.

Interpretation

The ENSO lowess showed a higher frequency pattern than the other two oceanic systems. The lowess I chose for the ENSO definitely smoothed over some higher frequencys spikes in the data, but overall there were about 5 year spans between the minima or maxima. The higher peaks in the ts plot can be explained by the periodogram of the ENSO which demonstrated a dominant frequency of one cycle per year, or an annual pattern. Wang et al. (2012) demonstrated that ENSO phases can last on the order of 3 to 5 years which corresponds to the second highest peak in the periodogram at frequency 0.2, or a period of 5 years. The lowess smoother for the PDO that looked the best to me highlighted a pattern of 30-50 years in between minima. This pairs well with the results of the periodogram showing dominant frequencies at 0.02 and 0.03 which correspond to periods of 30-50 years. There are higher frequencies at 0.1 and 0.2 which correspond to 3 and 5 year periods, though they are smaller than the frequency peak at 0.02. A recent study by Hu and Feng (2011) found that the circulation patterns set up by the PDO in the temperate latitudes could actually be affected by the strong anomalies of El Nino years. This may explain the higher frequency peaks in the PDO periodogram on the order of 3 to 5 years where the PDO is affected after El Nino years, creating a high frequency pattern with the multidecadal rhythym. Annual and mutlidecadal oceanic patterns in sea surface temperature anomalies change regional precipitation and circulation patterns all over the globe. The AMO lowess I chose showed a much lower frequency pattern, with peaks every 60-80 years. McCabe et al. 2012 found that the AMO showed circulation patterns of 60-80 years which is consistent with the period uncovered in the periodogram.The AMO periodogram had a dominant frequency at 0.01 through 0.03 that were by large greater than any other frequencies. This corresponds to a period of 50-100 years which is a much longer period than the other two patterns. There was also a smaller cycle present at higher frequencies of 0.12-0.14 which correspond to a period of 5-7 years. In addition, Hu and Feng (2011) found that the El Nino years strengthened the affect of AMO on North American Precipitation, which again can account for those higher frequency patterns in the multidecadal data. Hu and Feng then surmised that El Nino could be a dominant force in determining precipitation patterns due to the interannual cycle and its affect on the phases and strength of the PDO and AMO.

Carlos Duarte et al. 2012. Abrupt Climate Change in the Arctic. Nature Climate Change. Vol 2.

Hughes, M.K. and G. Funkhouse. 2003. Frequency-dependent climate signal in upper and lower forest border tree rings in the mountains of the Great Basin. Climatic Change 8:627-634.

LaMarche,V.C. 1974. Frequency-dependent relationships between tree-ring series along an ecological gradient and some dendroclimatic implications. Tree-Ring Bulletin 34:1-20.

LaMarch, V.C, D.A. Graybill, H.C. Fritts, and M.R. Rose.1984. Increasing atmospheric carbon dioxide: tree ring evidence for growth enhancement in natural vegetation. Science 225:219.

McCabe, G.J., T.R. Ault, B.I. Cook, J.L. Betancourt, and M.D. Schwartz.2012. Influences of the El Niño Southern Oscillation and the Pacific Decadal Oscillation on the timing of the North American spring. International Journal of Climatology 32:2301-2310.

Qi Hu and Song Feng. 2011.AMO- and ENSO-driven Summertime Circulation and Precipitation Variations in North America. Journal of the American Meteorlogical Society.

Schulze, E.D., H.A. Mooney, and E.L. Dunn. 1967. Wintertime photosynthesis of bristlecone pine.Journal of Ecology 48:1044-1047.

Tilling et al. 2015. Increased Arctic sea ice volume after anomalously low melting in 2013.Nature Geoscience 8:643-648.

Wang, B., J. Liu, H.J. Kim, P.J. Webster, S.Y. Yim, B. Xiang. 2012. Northern Hemisphere and summer monsoons intensified by mega ENSOs and AMO. PNAS 110:5347-5352