Assignment 2

Loading the required packages

library(sp)
library(raster) 
library(ncdf4)
library(rgdal)
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

Setting the working directory and file path

setwd("C:/MSc_CC/GAssignment2")
path <- file.path(getwd(),"Data")

Reading in the coastline Spatial Data Frame

world.coast <- readOGR(dsn = file.path(path,"ne_110m_coastline"), layer = "ne_110m_coastline" )
world.coast <- spTransform(world.coast, CRS( "+proj=longlat +datum=WGS84" ) )

Loading the Hadely SST data

nchad<-nc_open(file.path(path,'HadISST_sst.nc'))
print(nchad)

Extracting the latitude and longitude

lat<-ncvar_get(nchad, 'latitude')
lon<-ncvar_get(nchad, 'longitude') 

Extracting the time and formatting the units

timehad<-ncvar_get(nchad, 'time') 
tunitshad<-ncatt_get(nchad,"time",attname="units")
tustrhad<-strsplit(tunitshad$value, " ")
datehad<-as.character(as.Date(timehad,origin=unlist(tustrhad)[3]))

Extracting the Hadely SST data

ssthad <- ncvar_get(nchad, 'sst')

Replacing missing values with NAs

fillvaluehad <- ncatt_get(nchad,"sst","_FillValue")
ssthad[ssthad==fillvaluehad$value] <- NA
missvaluehad <- ncatt_get(nchad,"sst","missing_value")
ssthad[ssthad==missvaluehad$value] <- NA
ssthad[ssthad==-1000] <- NA

Creating an annual average of the Hadley SST data.

yearhad <- format(as.Date(datehad, format="%Y-%m-%d"),"%Y")
gmeanhad <- colMeans(ssthad, na.rm = TRUE, dims=2)
annmeanhad <- aggregate(gmeanhad,by=list(yearhad),FUN=mean,na.rm=TRUE)

Note: Annual average of Hadley sst (annmeanhad) runs from 1870-2017.

avssthad = rowMeans(ssthad,na.rm=TRUE,dims=2)
avssthad = rowMeans(ssthad,na.rm=FALSE,dims=2)

creating colours for plotting

colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)

levelplot(avssthad,col.regions = pal(100));

gridhad <- expand.grid(x=lon, y=lat)

gridhad$avssthad <- as.vector(avssthad)
levelplot(avssthad~x*y,gridhad,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SST'
) + 
  layer(sp.lines(world.coast))

Making annual averages from monthly data

yrshad <- annmeanhad$Group.1 
nyrhad <- length(yrshad)

assthad <- array(NA,c(dim(lon),dim(lat),nyrhad)) 

for (k in 1:nyrhad) {
  assthad[,,k] <- rowMeans(ssthad[,,yearhad==yrshad[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}

Adding annual averages to the previous data frame

gridhad$an_avssthad <- as.vector(rowMeans(assthad,na.rm=FALSE,dims=2))

Plotting this

levelplot(an_avssthad~x*y, data=gridhad,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Annually Averaged SST') + 
layer(sp.lines(world.coast)) 

Removing the global mean from each year

gmeanhad <- colMeans(assthad, na.rm = TRUE, dims=2)
for (k in 1:nyrhad){
  assthad[,,k]<-assthad[,,k]-matrix(gmeanhad[k],length(lon),length(lat))
}

Index 1: Atlantic Multidecadal Oscillation

Reading in the AMO data file

ncamo<-nc_open(file.path(path,'AMO.nc'))
print(ncamo)

Extracting the time and formatting the units

timeamo<-ncvar_get(ncamo, 'time') 
tunitsamo<-ncatt_get(ncamo,"time",attname="units")
tustramo<-strsplit(tunitsamo$value, " ")
dateamo<-as.character(as.Date(timeamo*365.25/12,origin=unlist(tustramo)[3]))

Extracting the AMO sst data

sstamo <- ncvar_get(ncamo, 'AMO')

Replacing the missing values with NAs

fillvalueamo <- ncatt_get(ncamo,"AMO","_FillValue")
sstamo[sstamo==fillvalueamo$value] <- NA
missvalueamo <- ncatt_get(ncamo,"AMO","missing_value")
sstamo[sstamo==missvalueamo$value] <- NA
sstamo[sstamo==3000000000000000000000000000000000] <- NA

creating an annual mean for the AMO index

amoyear <- format(as.Date(dateamo, format ="%Y-%m-%d"),"%Y")
amoannmean <- aggregate(sstamo,by=list(amoyear), FUN= mean, na.rm= TRUE)
head(amoannmean)

##   Group.1           x
## 1    1880  0.26576526
## 2    1881 -0.02418830
## 3    1882  0.06306856
## 4    1883  0.01907216
## 5    1884 -0.08529102
## 6    1885  0.04853270

colnames(amoannmean) <- c("Year", "AMO")

Note: AMO annual average (amoannmean) runs from 1880-2019.

TO BE WRITTEN ABOUT Renaming the column names in the annual average Hadley sst dataframe.

colnames(annmeanhad) <- c("Year", "Had")

Making the length of the annual averages of Hadley sst and AMO time series’ the same length.

library(dplyr)

## Warning: package 'dplyr' was built under R version 3.5.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

hadamomerge <- merge(annmeanhad, amoannmean, all=TRUE)
amohadmergedsubset <- hadamomerge %>% slice(11:148)

assthadsubset <- assthad[,,11:148]

Plotting this timeseries

plot(amohadmergedsubset$Year, amohadmergedsubset$AMO, type='l',xlab='Year',ylab='SST Anomaly',main='AMO')

Correlation of SSTA with AMV

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))


for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset[i,j,][!is.na(assthadsubset[i,j,])])>2){
      c.matrix[i,j] <- cor(assthadsubset[i,j,], as.numeric(amohadmergedsubset$AMO))
      p.vals <- cor.test(assthadsubset[i,j,], amohadmergedsubset$AMO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

Add to the same data frame as earlier for plotting

gridhad$corr <- as.vector(c.matrix)
gridhad$pval <- as.vector(t.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the correlation

levelplot(corr~x*y, data=gridhad , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with AMV')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

Linear regression of SSTA with AMO

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset[i,j,][!is.na(assthadsubset[i,j,])])>2){
      r.lm <- lm(assthadsubset[i,j,]~(as.numeric(amohadmergedsubset$AMO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

Adding this to the same data frame as earlier to plot

gridhad$reg <- as.vector(r.matrix)
gridhad$sig <- as.vector(s.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the linear regression for AMO and SSTA

levelplot(reg~x*y, data=gridhad , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SSTA with AMo', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Index 2: Pacific Decadal Oscillation

Reading in the PDO data

ncpdo<-nc_open(file.path(path,'PDO.nc'))
print(ncpdo)

Extracting the time and formatting the units

timepdo<-ncvar_get(ncpdo, 'time') 
tunitspdo<-ncatt_get(ncpdo,"time",attname="units")
tustrpdo<-strsplit(tunitspdo$value, " ")
datepdo<-as.character(as.Date(timepdo*365.25/12,origin=unlist(tustrpdo)[3]))

Extracting the PDO sst data

sstpdo <- ncvar_get(ncpdo, 'PDO')

Replacing the missing values with NAs

fillvaluepdo <- ncatt_get(ncpdo,"PDO","_FillValue")
sstpdo[sstpdo==fillvaluepdo$value] <- NA
missvaluepdo <- ncatt_get(ncpdo,"PDO","missing_value")
sstpdo[sstpdo==missvaluepdo$value] <- NA
sstpdo[sstpdo==3000000000000000000000000000000000] <- NA

Getting the annual mean of PDO

pdoyear <- format(as.Date(datepdo, format ="%Y-%m-%d"),"%Y")
pdoannmean <- aggregate(sstpdo,by=list(pdoyear), FUN= mean, na.rm= TRUE)
colnames(pdoannmean) <- c("Year", "PDO")

Note: the annual average of PDO (pdoannmean) runs from 1900-2017.

Merging the Hadley sst data and PDO to cut overlapping dates on each time series.

hadpdomerge <- merge(annmeanhad, pdoannmean, all=TRUE)

pdohadmergedsubset <- hadpdomerge %>% slice(31:148)

Subsetting the asst Hadley data to coincide with the PDO timeseries length.

assthadsubset2 <- assthad[,,31:148]

Plotting the PDO time series

plot(pdohadmergedsubset$Year, pdohadmergedsubset$PDO, type='l',xlab='Year',ylab='SST Anomaly',main='PDO')

Correlation of SSTA with PDO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset2[i,j,][!is.na(assthadsubset2[i,j,])])>2){
      c.matrix[i,j] <- cor(assthadsubset2[i,j,], as.numeric(pdohadmergedsubset$PDO))
      p.vals <- cor.test(assthadsubset2[i,j,], pdohadmergedsubset$PDO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

Adding this to the dataframe from earlier.

gridhad$corr <- as.vector(c.matrix)
gridhad$pval <- as.vector(t.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the Hadley sst with PDO correlation

levelplot(corr~x*y, data=gridhad , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with PDO')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

linear regression of SSTA with PDO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset2[i,j,][!is.na(assthadsubset2[i,j,])])>2){
      r.lm <- lm(assthadsubset2[i,j,]~(as.numeric(pdohadmergedsubset$PDO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

Adding this to the same data frame as earlier

gridhad$reg <- as.vector(r.matrix)
gridhad$sig <- as.vector(s.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plot the linear regression for Hadley sst and PDO

levelplot(reg~x*y, data=gridhad , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SSTA with ODI', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00005, col = "black"))

Index 3: El Nino Southern Oscillation

Reading in the ENSO data

ncenso<-nc_open(file.path(path,'ENSO.nc'))
print(ncenso)

Extracting the time and formatting the units

timeenso<-ncvar_get(ncenso, 'time')
tunitsenso<-ncatt_get(ncenso,"time",attname="units")
tustrenso<-strsplit(tunitsenso$value, " ")
dateenso<-as.character(as.Date(timeenso*365.25/12,origin=unlist(tustrenso)[3]))

Extracting the ENSO sst data

sstenso <- ncvar_get(ncenso, 'Nino4r')

Replacing the missing values with NAs

fillvalueenso <- ncatt_get(ncenso,"Nino4r","_FillValue")
sstenso[sstenso==fillvalueenso$value] <- NA
missvalueenso <- ncatt_get(ncenso,"Nino4r","missing_value")
sstenso[sstenso==missvalueenso$value] <- NA
sstenso[sstenso==3000000000000000000000000000000000] <- NA

Getting the annual mean of ENSO

ensoyear <- format(as.Date(dateenso, format ="%Y-%m-%d"),"%Y")
ensoannmean <- aggregate(sstenso,by=list(ensoyear), FUN= mean, na.rm= TRUE)
colnames(ensoannmean) <- c("Year", "ENSO")

Note: the annual avergae ENSO mean runs from 1854-2019.

Merging the Hadley sst data and ENSO to cut overlapping dates on each time series.

hadensomerge <- merge(annmeanhad, ensoannmean, all=TRUE)
ensohadmergedsubset <- hadensomerge %>% slice(17:164)

Subsetting the asst Hadley data to coincide with the ENSO timeseries length. The Hadley asst is already the same length as the ENSO data (from 1870-2017)/

assthadsubset3 <-  assthad

Plotting this time series

plot(ensohadmergedsubset$Year, ensohadmergedsubset$ENSO, type='l',xlab='Year',ylab='SST Anomaly',main='ENSO')

Correlation of SSTA with ENSO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset3[i,j,][!is.na(assthadsubset3[i,j,])])>2){
      c.matrix[i,j] <- cor(assthadsubset3[i,j,], as.numeric(ensohadmergedsubset$ENSO))
      p.vals <- cor.test(assthadsubset3[i,j,], ensohadmergedsubset$ENSO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

gridhad$corr <- as.vector(c.matrix)
gridhad$pval <- as.vector(t.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the correlation between Hadley sst and ENSO

levelplot(corr~x*y, data=gridhad , 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with ENSO')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

Linear regression of SSTA with ENSO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(assthadsubset3[i,j,][!is.na(assthadsubset3[i,j,])])>2){
      r.lm <- lm(assthadsubset3[i,j,]~(as.numeric(ensohadmergedsubset$ENSO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

gridhad$reg <- as.vector(r.matrix)
gridhad$sig <- as.vector(s.matrix)

sig <- subset(gridhad[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plot the linear regression for ENSO

levelplot(reg~x*y, data=gridhad , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SSTA with AMo', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Discussion of these patterns and discussion of the role that ocean circulation may have played in producing these patterns.

SST

AMO

In the North Atlantic where AMO operates, the correlation between SSTA and AMO is all generally higher than 0.3. The area highest correlation of SSTA with AMO is between the northern and central North Atlantic Ocean. This area of highest correlation, at around 0.7, extends from the west of Portugal, westwards to the centre and rises northwards towards the east coast of Canada and the north east coast of USA. Similarly, for the linear regression map of SSTA and AMO the highest regression value of 3 lies off the east coast of Canada and north eastern USA. Darker areas of high regression also then extend towards the centre to northern art of the Atlantic Ocean at values of around 2-3.

The Atlantic Multidecadal Oscillation (AMO) operates in the North Atlantic Ocean and can be recognised by phase changes in SST across the North Atlantic (McCarthy et al., 2015). These phase changes generally exist for around a decade. There is much debate in literature regarding the mechanisms which drive AMO phase changes (Trenary & Delsole, 2016). Some papers have pointed to the changes in SST are controlled by random atmospheric forcing, with no role of ocean dynamics; an argument based on using slab ocean models (Zhang, 2015). This paper found the key physical properties of AMO cannot be captures using just atmospheric mechanism in models, not in slab ocean models (Stanley, 2017). Ocean dynamic in models allow AMO to be reproduced (Delworth et al., 2007). Some papers argue ocean dynamics plays a dominant role in controlling NAO, particularly in the subpolar region where the highest correlations and regression is found in this study.

There may also be a link to variations in the subpolar gyre in causing AMO changes (Zhang et al., 2016). As the subpolar, and subtropical, gyres seen changes in circulation, the NAO and associated ocean circulation changes. This then has an impact on the SST in the North Atlantic and thus, the phase changes of AMO (McCarthy et al., 2015). There has been a link between AMOC and AMO. AMO a result of variation in AMOC (Delworth et al., 2007). Models have provided proof of this link by showing a clear multidecadal scale oscillation which links to AMOC (Mann et al., 2014). The AMO seems to be influenced by northward movement of heat by AMOC (Trenary & Delsole, 2016). A warm AMO phase coincides with a strengthened AMOC (wang & Zhang, 2013). The enhancing poleward heat transport= broadscale warming (Delworth et al., 2007). The AMOC is found to enhance AMO predictive skill in models- showing AMOC and AMO link. Most models capture multidecadal variability of SST and AMOC but to different amplitudes and frequencies (Trenary & Delsole, 2016).

PDO

The highest correlation of SSTA and PDO in the Pacific Ocean is all around the equator to the tropics of the Pacific Ocean, at around 0.8. Along the west coast of America also has high values of around 0.8. The west to centre of the North Pacific generally shows low correlation to SSTAs, with values of around 0.2- 0.4. The western to centre of the South Pacific also shows low correlation values of around 0.2-0.5. The linear regression of SSTA and PDO shows a similar pattern. The highest regression values are found across the entire basin along the equator, extending to the tropics north and south, with values of around 1. The lowest values are again found around the western to central parts of the North Pacific, an around the western to central parts of the South Pacific.

Many studies have established a link between PDO and the gyres in the North and South Pacific. Generally, a negative phase of PDO will coincide with warmer SST in the gyres, as well as warmer SSTs in the west of the Pacific, and colder in the east (Linsley et al., 2015). A warm phase of PDO is evident in the maps produced.A warm phase of PDO also coincides with a stronger and an expanded North Pacific subpolar gyre. This change in gyre shape is influences by a deeper Aleutian low. Studies also have shown SSTAs on a decadal timescale have a link to the gyre’s variability (Wills et al., date). However, studies have also shown in the Atlantic, when the gyre is contracted and smaller, there is less production of eddies and so less SSTAs. As such, a smaller gyre in the Pacific may also add to higher SSTAs as shown in the Atlantic (Qiu & Chen, 2005). Studies have also shown adjustments in ocean circulation which are linked to PDO phase changes. These readjustments are thought to be the cause of SSTAs which are associated with creating the PDO patterns shown in the maps (Wills et al., date).

Previously, the strength of the transport by the Kuroshio through the East China Sea has thought to be correlated to PDO variability (Andres et al., 2009). As seen in the Atlantic Ocean, the Gulfstream has aided producing SSTAs from a Northwards shift in its axis, meandering and a westward shift in the destabilisation point. This all aids in increases chance of circulation-ring interactions (Gawarkiewicz et al., 2018). The Kuroshio Extension may have a similar impact in adding to SSTAs associated with PDO. However, the link between the Kuroshio Extension and PDO has decreased since 1999, with changes in ocean circulation and gyres thought to be a driver of PDO variability (Wu et al., 2019). Finally, a link between PDO SSTAs in the Pacific and El Nino has been made (Linsley et al., 2015).

ENSO

The highest correlation values between SSTAs and ENSO can be seen all along the west coast of North and South America, with values of around 0.5-0.8. The highest values are along the equator at 0.8 and above. This high correlation, or SSTA warm tongue extends westward across the central Pacific, towards Papa New Guinea. The areas of lowest correlation of ENSO with SSTAs occur along the east coast of Australia, Pap New Guinea, and the Philippines, with values of around 0.3-0.6. Cooler, low correlation values occur in the North Pacific, extending from the east coast of Asia, towards the west coast of North America. The regression map shows a similar pattern with the higher regression values of SSTA and ENSO occurring along the west coast of the Americas, extending westwards towards the central Pacific. Values are generally between 0-2. The lower regression values occur along the east coast of Australia, Papa New Guinea and the Philippines, with negative values extending eastwards towards North America in the North Pacific.

The significance of these patterns from the correlation and linear regression suggests an El Nino year as warmer water is more towards the east and centre of the Pacific Ocean. This has this a significant impact on weather across the entire Pacific basin. While during normal years the east side of the Pacific has dry weather, the movement of such vast quantities of warm water east brings associated flooding, and heavy rainfall to South America (O’Hare et al., 2005). Impacts include reduced crop yields or flooded homes and businesses. On the western side, Australia, Papa New Guinea and the Philippines are subjected to reduced rainfall and drought conditions (Aguado & Burt, 2015). The upwelled cooler water resulting from the movement east of warm water, significantly reduces rainfall. ENSO events also have significant global impacts reaching further than the Pacific basin, such as influencing the Indian Monsoon and causing past ice storms across North America (Kumar et al., 2006) (Lecomte et al., 1998).

The drivers of ENSO, and then also the role of ocean circulation, is still an active area of research. While its understood naturally variability plays a key role, many have pointed to the role of climate change in driving events, both which can impact ocean circulation. A characteristic of ENSO is changes in the thermocline depth, which are thought to change under climate change. This may increase ENSO events, particularly central Pacific events (Turk et al., 2011). Circulation and upwelling variations may also influence ENSO events (Olderborgh, 2005).

Sea Level Pressure (slp)

Loading the SLP data

ncslp<-nc_open(file.path(path,'slp.mon.mean.nc'))

Extract and format the longitude and latitude

lat<-ncvar_get(ncslp, 'lat')
lon<-ncvar_get(ncslp, 'lon')
lon <- ifelse(lon>180, lon-360,lon) #changing the scale of lon values

Extracting the time and formating the units

timeslp<-ncvar_get(ncslp, 'time') 
tunitsslp<-ncatt_get(ncslp,"time",attname="units")
tustrslp<-strsplit(tunitsslp$value, " ")
dateslp<-as.character(as.Date(timeslp/24,origin=unlist(tustrslp)[3]))

note:ssthad=slp

Extracting the slp data

slp <- ncvar_get(ncslp, 'slp')

Replacing missing values with NAs

fillvalueslp <- ncatt_get(ncslp,"slp","_FillValue")
slp[slp==fillvalueslp$value] <- NA
missvalueslp <- ncatt_get(ncslp,"slp","missing_value")
slp[slp==missvalueslp$value] <- NA
slp[slp==-1000] <- NA

Creating annual averages

yearslp <- format(as.Date(dateslp, format="%Y-%m-%d"),"%Y")
gmeanslp <- colMeans(slp, na.rm = TRUE, dims=2)
annmeanslp <- aggregate(gmeanslp,by=list(yearslp),FUN=mean,na.rm=TRUE)

avsstslp = rowMeans(slp,na.rm=TRUE,dims=2)
avsstslp = rowMeans(slp,na.rm=FALSE,dims=2)

levelplot(avsstslp,col.regions = pal(100));

gridslp <- expand.grid(x=lon, y=lat)

gridslp$avsstslp <- as.vector(avsstslp)

Plotting the SLP

levelplot(avsstslp~x*y,gridslp,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SLP'
) + 
  layer(sp.lines(world.coast))

Making annual averages from monthly data

yrsslp <- annmeanslp$Group.1
nyrslp <- length(yrsslp)

asstslp <- array(NA, c(dim(lon),dim(lat),nyrslp)) 

for (k in 1:nyrslp) {
  asstslp[,,k] <- rowMeans(slp[,,yearslp==yrsslp[k]],na.rm=FALSE,dims=2)
}

Removing the global mean from each year

gridslp$an_avsstslp <- as.vector(rowMeans(asstslp,na.rm=FALSE,dims=2))

Plotting this

levelplot(avsstslp~x*y, data=gridslp,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') + 
  layer(sp.lines(world.coast)) 

Removing the global mean from each year

gmeanslp <- colMeans(asstslp, na.rm = TRUE, dims=2)
for (k in 1:nyrslp){
  asstslp[,,k]<-asstslp[,,k]-matrix(gmeanslp[k],length(lon),length(lat))
}

Index 1: Atlantic Multidecadal Oscillation

Reading in the AMO data file

ncamo<-nc_open(file.path(path,'AMO.nc'))
print(ncamo)

Extracting the time and formatting the units

timeamo<-ncvar_get(ncamo, 'time') 
tunitsamo<-ncatt_get(ncamo,"time",attname="units")
tustramo<-strsplit(tunitsamo$value, " ")
dateamo<-as.character(as.Date(timeamo*365.25/12,origin=unlist(tustramo)[3]))

Extracting the AMO sst data

sstamo <- ncvar_get(ncamo, 'AMO')

Replacing the missing values with NAs

fillvalueamo <- ncatt_get(ncamo,"AMO","_FillValue")
sstamo[sstamo==fillvalueamo$value] <- NA
missvalueamo <- ncatt_get(ncamo,"AMO","missing_value")
sstamo[sstamo==missvalueamo$value] <- NA
sstamo[sstamo==3000000000000000000000000000000000] <- NA

creating an annual mean for the AMO index

amoyear <- format(as.Date(dateamo, format ="%Y-%m-%d"),"%Y")
amoannmean <- aggregate(sstamo,by=list(amoyear), FUN= mean, na.rm= TRUE)
head(amoannmean)

##   Group.1           x
## 1    1880  0.26576526
## 2    1881 -0.02418830
## 3    1882  0.06306856
## 4    1883  0.01907216
## 5    1884 -0.08529102
## 6    1885  0.04853270

colnames(amoannmean) <- c("Year", "AMO")

Note: AMO annual average (amoannmean) runs from 1880-2019.

Renaming the column names in the annual average SLP dataframe in order to merge

colnames(annmeanslp) <- c("Year", "SLP")

Making the length of the annual averages of Hadley sst and AMO time series’ the same length.

library(dplyr)
slpamomerge <- merge(annmeanslp, amoannmean, all=TRUE)
amoslpmergedsubset <- slpamomerge %>% slice(69:139)

asstslpsubset <- slp[,,69:139]

Plotting this timeseries

plot(amoslpmergedsubset$Year, amoslpmergedsubset$AMO, type='l',xlab='Year',ylab='SLP Anomaly',main='AMO')

Correlation of SLP with AMO

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))


for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset[i,j,][!is.na(asstslpsubset[i,j,])])>2){
      c.matrix[i,j] <- cor(asstslpsubset[i,j,], as.numeric(amoslpmergedsubset$AMO))
      p.vals <- cor.test(asstslpsubset[i,j,], amoslpmergedsubset$AMO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

Add to the same data frame as earlier for plotting

gridslp$corr <- as.vector(c.matrix)
gridslp$pval <- as.vector(t.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the correlation

levelplot(corr~x*y, data=gridslp , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with AMO')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

Linear regression of SLP with AMO

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset[i,j,][!is.na(asstslpsubset[i,j,])])>2){
      r.lm <- lm(asstslpsubset[i,j,]~(as.numeric(amoslpmergedsubset$AMO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

Adding this to the same data frame as earlier to plot

gridslp$reg <- as.vector(r.matrix)
gridslp$sig <- as.vector(s.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the linear regression for AMO and SLP

levelplot(reg~x*y, data=gridslp , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SLP with AMO', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Index 2: Pacific Decadal Oscillation

Reading in the PDO data

ncpdo<-nc_open(file.path(path,'PDO.nc'))
print(ncpdo)

Extracting the time and formatting the units

timepdo<-ncvar_get(ncpdo, 'time') 
tunitspdo<-ncatt_get(ncpdo,"time",attname="units")
tustrpdo<-strsplit(tunitspdo$value, " ")
datepdo<-as.character(as.Date(timepdo*365.25/12,origin=unlist(tustrpdo)[3]))

Extracting the PDO sst data

sstpdo <- ncvar_get(ncpdo, 'PDO')

Replacing the missing values with NAs

fillvaluepdo <- ncatt_get(ncpdo,"PDO","_FillValue")
sstpdo[sstpdo==fillvaluepdo$value] <- NA
missvaluepdo <- ncatt_get(ncpdo,"PDO","missing_value")
sstpdo[sstpdo==missvaluepdo$value] <- NA
sstpdo[sstpdo==3000000000000000000000000000000000] <- NA

Getting the annual mean of PDO

pdoyear <- format(as.Date(datepdo, format ="%Y-%m-%d"),"%Y")
pdoannmean <- aggregate(sstpdo,by=list(pdoyear), FUN= mean, na.rm= TRUE)
colnames(pdoannmean) <- c("Year", "PDO")

Note: the annual average of PDO (pdoannmean) runs from 1900-2017.

Merging the SLP data and PDO to cut overlapping dates on each time series.

slppdomerge <- merge(annmeanslp, pdoannmean, all=TRUE)

pdoslpmergedsubset <- slppdomerge %>% slice(31:148)
pdoslpmergedsubset <- slppdomerge %>% slice(49:118)

Subsetting the SLP data to coincide with the PDO timeseries length.

asstslpsubset2 <- slp[,,1:70]

Plotting the PDO time series

plot(pdoslpmergedsubset$Year, pdoslpmergedsubset$PDO, type='l',xlab='Year',ylab='SLP Anomaly',main='PDO')

Correlation of SLP with PDO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset2[i,j,][!is.na(asstslpsubset2[i,j,])])>2){
      c.matrix[i,j] <- cor(asstslpsubset2[i,j,], as.numeric(pdoslpmergedsubset$PDO))
      p.vals <- cor.test(asstslpsubset2[i,j,], pdoslpmergedsubset$PDO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

Adding this to the dataframe from earlier.

gridslp$corr <- as.vector(c.matrix)
gridslp$pval <- as.vector(t.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the SLP with PDO correlation

levelplot(corr~x*y, data=gridslp , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with PDO')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

linear regression of SLP with PDO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset2[i,j,][!is.na(asstslpsubset2[i,j,])])>2){
      r.lm <- lm(asstslpsubset2[i,j,]~(as.numeric(pdoslpmergedsubset$PDO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

Adding this to the same data frame as earlier

gridslp$reg <- as.vector(r.matrix)
gridslp$sig <- as.vector(s.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plot the linear regression for Hadley sst and PDO

levelplot(reg~x*y, data=gridslp , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SLP with PDO', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00005, col = "black"))

Index 3: El Nino Southern Oscillation

Reading in the ENSO data

ncenso<-nc_open(file.path(path,'ENSO.nc'))
print(ncenso)

Extracting the time and formatting the units

timeenso<-ncvar_get(ncenso, 'time')
tunitsenso<-ncatt_get(ncenso,"time",attname="units")
tustrenso<-strsplit(tunitsenso$value, " ")
dateenso<-as.character(as.Date(timeenso*365.25/12,origin=unlist(tustrenso)[3]))

Extracting the ENSO sst data

sstenso <- ncvar_get(ncenso, 'Nino4r')

Replacing the missing values with NAs

fillvalueenso <- ncatt_get(ncenso,"Nino4r","_FillValue")
sstenso[sstenso==fillvalueenso$value] <- NA
missvalueenso <- ncatt_get(ncenso,"Nino4r","missing_value")
sstenso[sstenso==missvalueenso$value] <- NA
sstenso[sstenso==3000000000000000000000000000000000] <- NA

Getting the annual mean of ENSO

ensoyear <- format(as.Date(dateenso, format ="%Y-%m-%d"),"%Y")
ensoannmean <- aggregate(sstenso,by=list(ensoyear), FUN= mean, na.rm= TRUE)
colnames(ensoannmean) <- c("Year", "ENSO")

Note: the annual avergae ENSO mean runs from 1854-2019.

Merging the SLP data and ENSO to cut overlapping dates on each time series.

slpensomerge <- merge(annmeanslp, ensoannmean, all=TRUE)
ensoslpmergedsubset <- slpensomerge %>% slice(95:165)

Subsetting the SLP data to coincide with the ENSO timeseries length.

asstslpsubset3 <-  asstslp

Plotting this time series

plot(ensoslpmergedsubset$Year, ensoslpmergedsubset$ENSO, type='l',xlab='Year',ylab='SLP Anomaly',main='ENSO')

Correlation of SLP with ENSO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset3[i,j,][!is.na(asstslpsubset3[i,j,])])>2){
      c.matrix[i,j] <- cor(asstslpsubset3[i,j,], as.numeric(ensoslpmergedsubset$ENSO))
      p.vals <- cor.test(asstslpsubset3[i,j,], ensoslpmergedsubset$ENSO)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

gridslp$corr <- as.vector(c.matrix)
gridslp$pval <- as.vector(t.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plotting the correlation between SLP and ENSO

levelplot(corr~x*y, data=gridslp , 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with ENSO')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.00001, col = "black"))

Linear regression of SLP with ENSO

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asstslpsubset3[i,j,][!is.na(asstslpsubset3[i,j,])])>2){
      r.lm <- lm(asstslpsubset3[i,j,]~(as.numeric(ensoslpmergedsubset$ENSO)))
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

gridslp$reg <- as.vector(r.matrix)
gridslp$sig <- as.vector(s.matrix)

sig <- subset(gridslp[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

Plot the linear regression for ENSO

levelplot(reg~x*y, data=gridslp , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste('Regression of SLP with AMO', lon, ',Lat=', lat)) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Discussion of these patterns and discussion of the role that ocean circulation may have played in producing these patterns.

SLP

AMO

The correlation and linear regression map both show a similar pattern. Focusing on the Atlantic Ocean as AMO operates here, the highest correlation and regression values tends to be in the central North Atlantic. Particularly off the west coast of Africa and Portugal, correlation values are around 0.2 and regression values are around 3. The lower correlation and regression values occur in the northern most parts of the North Atlantic and the southern most parts of the North Atlantic. Correlation values tend to be around 0.1 and regression values of around -1.

The phase changes of AMO can have significant impact on the global weather. For example, a lower mean sea level pressure, which is linked to the positive AMO phase has significant implications for precipitation. This phase and associated SLP causes a reduction in precipitation over North America, while Europe will have an increase of precipitation, particularly in the North and West (Knight et al., 2006).

PDO

The highest correlation and regression values are found in the eastern to central parts of the Pacific basin. Correlation values here are around 0.2-0.3, while regression values are around 0- 1. The lowest correlations and regression values of SLP to PDO is found in the north eastern region, off the west coast of North America. Values of correlation and regression are low also to the west of the Pacific basin, or off the east coast of Asia. Correlation values are around -0.3 and regression values are around -1.

The sea level pressure pattern is a key characteristic to identifying PDO. PDO and El Nino show a similar SLP pattern. During positive phases of PDO, the western side of the Pacific basin shows higher SLP. The patterns on the map show this well. The eastern side shows a lower SLP which stretches across towards the central Pacific (Yeh et al., 2018). Potentially, El Nino events and associated drivers of events in the ocean circulatory system may influence PDO.

ENSO

SLP and ENSO correlation and regression values is highest to the west of the Pacific basin. Correlation values are around 0.2-0.4 and regression values are around 0-1. The majority of the Pacific basin appears to have a negative correlation and regression value. From the centre to east of the basin correlation values are around 0.4 - 0.6. Regression values in this area are around 0-2.

The patterns on the maps are to be expected. During an El Nino event, sea level pressure becomes lower over the eastern and central regions of the Pacific Ocean basin than what is typical during a non-El Nino year (Yeh et al., 2018). Whereas on the western side, higher pressure operates during al El Nino year which the maps capture well (O’Hare et al., 2005). The map therefor shows clearly the presence of EL Nino events based on SLP. This has the associated impact of then weakening the trade winds (Yeh et al., 2018). Weaker trade winds are a key cause of El Nino events, showing SLP can influence the occurrence of events.

References

Aguado, E. & Burt, J. E. (2015) Understanding Weather and Climate. Pearson Education Limited. Andres, M., Park, J., Wimbush, M., Zhu, X., Nakamura, H., Kim, K. & Chang, K. (2009) “Manifestation of the Pacific Decadal Oscillation in the Kuroshio” Geophysical Research Letters. Vol 36. Accessed 28/04/19 at https://www.whoi.edu/cms/files/Andres_et_al_2009_GRL_102086.pdf.

Delworth, T. L., Zeng, F., Zhang, L., Zhang, R., Vecchi, G. A. & Yang, X. (2017) “The Central Role of Ocean Dynamics in Connecting the North Atlantic Oscillation to the Extratropical Component of the Atlantic Multidecadal Oscillation” Journal of Climate. Accessed 29/04/19 at https://journals.ametsoc.org/doi/full/10.1175/JCLI-D-16-0358.1

Gawarkiewicz, G., Todd, R. E., Zhang, W., Partida, J., Gangopadhyay, A., Monim, M., Fratantoni, P., Mercer, A. M. & Dent, M. (2018) “The Changing Nature of Shelf-Break Exchange Revealed by the OOI Pioneer Array” Oceanography. Vol 31.

Knight, J. R., Folland, C. K. & Scaife, A. A. (2006) “Climate Impacts of the Atlantic Multidecadal Oscillation” Geophysical research letters. Vol 33.

Kumar, K. K., Rajagopalam, B., Hoerling, M., Bates, G. & Cane, M. (2006) “Unravelling the Mystery of the Indian Monsoon Failure During El Niño” Science., Vol 314, Pg. 115-118. Lecomte, E. L., Pang, A. W. & Russell, J. W. (1998) “Ice Storm ’98” ICLP Research Paper Series- No.1.

Linsley, B. K., Wu. H. C., Dassie, E. P. & Schrag, D. P. (2015) “Decadal changes in South Pacific sea surface temperatures and the relationship to the Pacific decadal oscillation and upper ocean heat content” Geophysical Research Letters. Vol 42. Accessed 27/04/19 at https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015GL063045.

Mann, M. E., Steinman, B. A. & Miller, S. K. (2014) “On forced temperature changes, internal variability, and the AMO” Geophysical Research Letters. Vol 41. Accessed 27/04/19 at https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014GL059233

McCarthy, G., Haigh, I. D., Hirschi, J. J. M., Grist, J. P. & Smeed, D. A. (2015) “Ocean impact on decadal Atlantic climate variability revealed by sea-level observations” Nature. Vol 521. Accessed 28/04/19 at https://www.nature.com/articles/nature14491

O’Hare, G., Sweeney, J. & Wilby, R. (2005) Weather, Climate and Climate Change - Human Perspectives. Pearson Education Limited . Oldenborgh, G. J. V., Philip, S. Y. & Collins, M. (2005) “El Niño in a changing climate: a multi-model study” Ocean Science Discussions., Vol 2, Pg.269-280. Qiu, B. & Chen, S. (2005) “Variability of the Kuroshio Extension Jet, Recirculation Gyre and Mesoscale Eddies on Decadal Time Scales” Journal of Physical Oceanography. Vol 35.

Stanley, S. (2017) “Ocean dynamics may drive North Atlantic temperature anomalies” Eos. Vol 98. Accessed 27/04/19 at https://eos.org/research-spotlights/ocean-dynamics-may-drive-north-atlantic-temperature-anomalies

Trenary, L. & Delsole, T. (2016) “Does the Atlantic Multidecadal Oscillation get its Predictability from the Atlantic Meridional overturning Circulation?” Journal of Climate. Accessed 28/04/19 at https://journals.ametsoc.org/doi/full/10.1175/JCLI-D-16-0030.1 Turk, D., Meinen, C. S., Antoiine, D., McPhaden, M. J. & Lewis, M. R. (2011) “Implications of Changing El Niño patterns for Biological Dynamics in the Equatorial Pacific Ocean” Geophysical Research Letters., Vol 38, Pg.1-6. Wang, C. & Zhang, L. (2013) “Multidecadal Ocean Temperature and Salinity Variability in the Tropical North Atlantic: Linking with the AMO, AMOC, and Subtropical Cell” Journal of Climate. Accessed 27/04 at https://journals.ametsoc.org/doi/full/10.1175/JCLI-D-12-00721.1.

Wills. R. C. J., Battisti, D. S., Proistosescu, C., Thompson, L., Hartmann, D. L. & Armour, K. C. (2019) “Ocean Circulation Signatures of North Paci???c Decadal Variability” Geophysical Research Letters. Accessed 28/04/19 at https://atmos.washington.edu/~dennis/2018_PDO.pdf.

Wu, C., Wang, Y. & Chao, S. (2019) “Disassociation of the Kuroshio Current with the Paci???c Decadal Oscillation Since 1999” Remote Sensing. Vol 11. Accessed 26/04/19 at https://www.mdpi.com/2072-4292/11/3/276?type=check_update&version=1

Yeh, S., Cai, W., Min, S., McPhaden, M. J., Dommenget, D., Dewitte, B., Collins, M., Ashok, K., An, S., Yim, B. & Kug, J. (2018) “ENSO Atmospheric Teleconnections and Their Response to Greenhouse Gas Forcing” Sea State and boundary Layer Physics of the Emerging Artic Ocean. Vol 56. Accessed 29/04/19 at https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017RG000568.

Zhang, R., Sutton, R., Danabasoglu, G., Delworth, T., Kim, W. M., Robson, J. & Yeager, S. G. (2016) “Commenton”The Atlantic Multidecadal Oscillationwithout a role for ocean circulation" Science. Vol 352. Accessed 29/04/19 at https://science.sciencemag.org/content/sci/352/6293/1527.1.full.pdf

Zhang, R. (2015) “On the persistence and coherence of subpolar sea surface temperature and salinity anomalies associated with the Atlantic multidecadal variability” Geophysical Research letters. Vol 44. Accessed 27/04/19 at https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/2017GL074342?purchase_referrer=agupubs.onlinelibrary.wiley.com&tracking_action=preview_click&r3_referer=wol&show_checkout=1.