Influence of global Warming on Western North Pacific tropical cyclone intensities during 2015

Se-hwan Yang, Nam-Young Kang, James B. Elsner and Young-Shin Chun

Abstract

The climate of 2015 was characterized by a strong El Nino, global warmth, and record-setting tropical cyclone (TC) intensity for western North Pacific typhoons. In this study, the highest TC intensity in 32 years (1984–2015) is shown to be a consequence of above normal TC activity—following natural internal variation—and greater efficiency of intensity. The efficiency of intensity (EINT) is termed the `blasting’ effect and refers to typhoon intensification at the expense of occurrence. Statistical models show that the EINT is mostly due to the anomalous warmth in the environment indicated by global mean sea-surface temperature. In comparison, the EINT due to El Ni~no is negligible. This implies that the record-setting intensity of 2015 might not occurred without environmental warming and suggests that a year with even greater TC intensity is possible in the near future when above normal activity coincides with another record EINT due to continued multidecadal warming.

Figure 1. Time series of NSOI & GMSST

##
rm(list=ls())

par(mfrow=c(1,1),mar=c(5,5,5,6))

#------------set option & read data---------------------------
year=1951:2015

NSOI=-1*read.table('Data/SOIann_jjason_1951to2015.dat',T)[(year[1]-1950):(rev(year)[1]-1950),1]  
GMSST=read.table('Data/GMSSTann_jjason_1854to2015.dat',T)[(year[1]-1853):(rev(year)[1]-1853),1]

#-------------plot--------------------------------------------
par(mar=c(4, 4, 4, 5))
plot(-999, -999, xlim=c(1950, 2016), ylim=c(-4, 4), axes=FALSE, xlab='', ylab='')
abline(h=0, lty=3)

polygon(c(1984,2015,2015,1984),c(10,10,-10,-10),border=NA,col='light yellow')
lines(year,scale(NSOI),col='purple',lwd=3) 
lines(year,scale(GMSST),col='orange',lwd=3) 

axis(1, at = seq(1950, 2020, 10), labels=seq(1950, 2020, 10), cex.axis=1)
axis(2, at = seq(-5,5,1), labels=seq(-5,5,1), las=1, cex.axis=1, col=1, col.axis=1)
gmsstlev=seq(17,19,0.2) ; gmsstsdlev=(gmsstlev-mean(GMSST))/sd(GMSST)
gmsstlev2=c('17.0','17.2','17.4','17.6','17.8','18.0','18.2','18.4','18.6','19.8','19.0')
axis(4, at = gmsstsdlev, labels=gmsstlev2, las=1, cex.axis=1, col=1, col.axis=1)

mtext('Year', side=1, line=2.6, cex=1)
mtext('Standardized negative SOI (s.d.)', side=2, line=2.5, cex=.95, col=1)
text(par("usr")[2] +10, 2.5, srt=-90, adj = 0, labels = expression(Global~mean~SST~(degree~C))
      , cex=.95, col=1, xpd = TRUE) 

legend('topleft',c('negative SOI','Global mean SST', 'scope of the investigation'),bty='n',
pch=c(-1,-1,15), pt.cex=1.5, col=c('purple','orange','light yellow'),lwd=3,cex=0.8)
box()

Fig. 1. Annual variation of the standardized values of global mean SST and negative SOI over the 65 years (1951–2015). JJASON observations are averaged and represented as annual values to represent the boreal summer environment for active TCs. Shaded period is the scope for this study.

Figure 2

require(png)

## Loading required package: png

## Warning: package 'png' was built under R version 3.3.2

require(grid)

## Loading required package: grid

img <- readPNG('diagram4Fig2.png')
grid.raster(img)

Fig. 2. The schematic of TC climate framework. The diagonal line labeled ACT indicates the in-phase relationship between FRQ and INT. The diagonal line labeled EINT indicates the out-of-phase relationship between FRQ and INT. Modified from Kang and Elsner (2012a).

Figure 3. Hovmoller diagram of TC climate

rm(list=ls())

# Read data
#------------------------------------------------------------------------------
INTjt=read.table('Data/intjtwc_jjason.txt',T)[,1]
INTjm=read.table('Data/intjma_jjason.txt',T)[,1]
FRQjt=read.table('Data/frqjtwc_jjason.txt',T)[,1]
FRQjm=read.table('Data/frqjma_jjason.txt',T)[,1]

#composite
INTcomp=((scale(INTjt)+scale(INTjm))/sqrt(2))
FRQcomp=((scale(FRQjt)+scale(FRQjm))/sqrt(2))

#-------------loop------------------------
library(akima) ; library(gtools)  

## Warning: package 'akima' was built under R version 3.3.2

## Warning: package 'gtools' was built under R version 3.3.2

layout(cbind(1,2))
par(mar=c(5,2,5,2),cex.axis=0.9)

vx=cbind(INTcomp,INTjt,INTjm)
vy=cbind(FRQcomp,FRQjt,FRQjm) 

for (vn in 1:3)
{
 # Scan 
  theta=seq(0,360,by=22.5)
  var.pass=NULL
  for (i in theta)
  {VR=cos(i*pi/180)*scale(vx[,vn])+sin(i*pi/180)*scale(vy[,vn])
    Fn=ecdf(VR) ; VR=Fn(VR)
    var.pass=rbind(var.pass,VR)
   }
  if (vn==1) var.com=var.pass
  if (vn==2) var.jtwc=var.pass
  if (vn==3) var.jma=var.pass
}

# plot result
  # ------------------------------------------------------------
  A1=rep(theta,each=32)  # 32 years (1984-2015)
  A2=rep(1:32,length(theta))       # 32 years (1984-2015) 
  B1=NULL ;  for(i in 1:length(theta)) {B1=c(B1,as.vector(var.com[i,]))} ; A3=B1  # for merged best track data
  B2=NULL ;  for(i in 1:length(theta)) {B2=c(B2,as.vector(var.jtwc[i,]))} ; A4=B2 # for JTWC best track data 
  B3=NULL ;  for(i in 1:length(theta)) {B3=c(B3,as.vector(var.jma[i,]))} ; A5=B3 # for JMA best track data 

par(mar=c(5,5,5,2),cex.axis=1.2)
  image(interp(A1,A2,A3, xo=seq(min(A1),max(A1),length=500),yo=seq(min(A2),max(A2),length=500),duplicate='mean')
        ,xlim=c(0,360),ylim=c(32,1),axes=F,xlab='',ylab='',col=c('#FFFFCC','#FFFF99','#FFFF66','#FFFF33','#FFCC00','#FF9900','#FF6600','#FF3300','#FF0000'))
  
  contour(interp(A1,A2,A3,xo=seq(min(A1),max(A1),length=100),yo=seq(min(A2),max(A2),length=100),duplicate='mean')        ,lwd=0.5, add = TRUE, col=colors()[280], labcex=1,levels=c(0,0.2,0.4,0.6,0.8))
  axis(1, at = c(0,45,90,135,180,225,270,315,360),labels= c('I','A','F','-E','-I','-A','-F','E','I'))
  axis(3, at = c(0,45,90,135,180,225,270,315,360),labels= c(0,'',90,'',180,'',270,'',360))
  axis(2, at = c(2,7,12,17,22,27,32),labels= c(1985,1990,1995,2000,2005,2010,2015),las=1)
  axis(4, at = c(2,7,12,17,22,27,32),labels= c('','','','','','',''),las=1)
  
  mtext('Direction',side=1,line=2.8, cex=1.5)
  mtext('Angle',side=3,line=2.8, cex=1.5)
  mtext('Year',side=2,line=3.4, cex=1.5)
  text(par("usr")[2] +120, 18, srt=-90, adj = 0, labels = expression(paste('  ')), cex=0.9, xpd = TRUE) 
  mtext('(a)',side=3,line=3,cex=1.9,adj=-0.17)
  box()

   #----------------------------------
   image(interp(A1,A2,A4-A5, xo=seq(min(A1),max(A1),length=500),yo=seq(min(A2),max(A2),length=500),duplicate='mean')
        ,xlim=c(0,360),ylim=c(32,1),axes=F,xlab='',ylab='' ,col=c('#00CCFF','#33FFFF','#66FFFF','#99FFFF','#CCFFFF','#F0FFFF','#FFFFF0','#FFFFCC','#FFFF99','#FFFF66','#FFFF33','#FFCC00'))
  
  contour(interp(A1,A2,A4-A5
                 ,xo=seq(min(A1),max(A1),length=100),yo=seq(min(A2),max(A2),length=100),duplicate='mean')       
          ,lwd=0.5, add = TRUE, col=colors()[280], labcex=1,levels=c(-0.4,-0.2,0,0.2,0.4))
  axis(1, at = c(0,45,90,135,180,225,270,315,360),labels= c('I','A','F','-E','-I','-A','-F','E','I'))
  axis(3, at = c(0,45,90,135,180,225,270,315,360),labels= c(0,'',90,'',180,'',270,'',360))
  axis(2, at = c(2,7,12,17,22,27,32),labels= c(1985,1990,1995,2000,2005,2010,2015),las=1)
  axis(4, at = c(2,7,12,17,22,27,32),labels= c('','','','','','',''),las=1)
  
  mtext('Direction',side=1,line=2.8, cex=1.5)
  mtext('Angle',side=3,line=2.8, cex=1.5)
  text(par("usr")[2] +120, 18, srt=-90, adj = 0, labels = expression(paste('  ')), cex=0.9, xpd = TRUE) 
  mtext('(b)',side=3,line=3,cex=1.9,adj=-0.17)
  box()

Fig. 3. Hovmoller diagram of (a) the annual variations of ranked probabilities in the merged WNP TC climate indicators during JJASON from JTWC and JMA best-track data (1984–2015), and (b) the difference between the ranked probabilities from the two observations. Probabilities are computed around the phase of the plane by the variable and principal component axes indicating INT, ACT, FRQ and EINT. INT, ACT, FRQ and EINT are denoted as I, A, F and E. Angles of 0, 45, 90, 135, 180, 225, 270, 315, and 360 degrees indicate I, A, F, -E, -I, -A, -F, E and I, respectively. Contours denote the ranked probability, which represents the probability level from an empirical cumulative density of the annual values. Climate indicators (lower abscissa) and equivalent angles (upper abscissa) are shown along the horizontal axis. Modified from Kang and Elsner (2012b).

Table 1. Trend of WNP TC activity

rm(list=ls())

#------------read data-------------------------------------
year=1984:2015

PDIjm=read.table('Data/pdijma_jjason.txt',T)[,1]
PDIjt=read.table('Data/pdijtwc_jjason.txt',T)[,1]

ACEjm=read.table('Data/acejma_jjason.txt',T)[,1]
ACEjt=read.table('Data/acejtwc_jjason.txt',T)[,1]

INTjt=read.table('Data/intjtwc_jjason.txt',T)[,1]*0.51444  #kt to m/s
INTjm=read.table('Data/intjma_jjason.txt',T)[,1]*0.51444  #kt to m/s
FRQjt=read.table('Data/frqjtwc_jjason.txt',T)[,1]
FRQjm=read.table('Data/frqjma_jjason.txt',T)[,1]

#composite
PDIcomp=(scale(PDIjt)+scale(PDIjm))/sqrt(2)
ACEcomp=(scale(ACEjt)+scale(ACEjm))/sqrt(2)

INTcomp=((scale(INTjt)+scale(INTjm))/sqrt(2))
FRQcomp=((scale(FRQjt)+scale(FRQjm))/sqrt(2))
ACTcomp=(scale(INTcomp)+scale(FRQcomp))/sqrt(2)

A=summary(lm(scale(PDIcomp)~year))$coefficients[2,1:4]  #PDI trend
B=summary(lm(scale(ACEcomp)~year))$coefficients[2,1:4]  #ACE trend
C=summary(lm(scale(ACTcomp)~year))$coefficients[2,1:4]  #ACT trend

PDIres=round(A,3)
ACEres=round(B,3)
ACTres=round(C,3)

data.frame(PDIres,ACEres,ACTres)

##            PDIres ACEres ACTres
## Estimate   -0.018 -0.027 -0.026
## Std. Error  0.019  0.019  0.019
## t value    -0.959 -1.457 -1.357
## Pr(>|t|)    0.345  0.156  0.185

#                  PDIres ACEres ACTres
# Estimate   -0.018 -0.027 -0.026
# Std. Error  0.019  0.019  0.019
# t value    -0.959 -1.457 -1.357
# Pr(>|t|)    0.345  0.156  0.185

Tab. 1. Trend of WNP TC activity during JJASON over the 32 years (1984–2015). PDI, ACE and ACT represent the merged variations from JTWC and JMA best-track data. P-val is computed under the null hypothesis of no trend.

Figure 4. Contribution of EINT to 2015 INT

rm(list=ls())

colset=c('black','#ff7300','blue','green','purple','#e1ff00', '#ffb400', '#9933cc30', '#5500ff')

#------------read data-------------------------------------
year=1984:2015

NSOI=-1*read.table('Data/SOIann_jjason_1951to2015.dat',T)[(year[1]-1950):(rev(year)[1]-1950),1]  
GMSST=read.table('Data/GMSSTann_jjason_1854to2015.dat',T)[(year[1]-1853):(rev(year)[1]-1853),1]
sNSOI=scale(NSOI)[,1]
sGMSST=scale(GMSST)[,1]

INTjt=read.table('Data/intjtwc_jjason.txt',T)[,1]*0.51444  #kt to m/s
INTjm=read.table('Data/intjma_jjason.txt',T)[,1]*0.51444  #kt to m/s
FRQjt=read.table('Data/frqjtwc_jjason.txt',T)[,1]
FRQjm=read.table('Data/frqjma_jjason.txt',T)[,1]
ACEjt=read.table('Data/acejtwc_jjason.txt',T)[,1]
ACEjm=read.table('Data/acejma_jjason.txt',T)[,1]

#composite
INTcomp=((scale(INTjt)+scale(INTjm))/sqrt(2))
FRQcomp=((scale(FRQjt)+scale(FRQjm))/sqrt(2))
ACTcomp=(scale(INTcomp)+scale(FRQcomp))/sqrt(2)
EINTcomp=(scale(INTcomp)-scale(FRQcomp))/sqrt(2)
ACEcomp=(scale(ACEjt)+scale(ACEjm))/sqrt(2)

#------------plot  role of ACT & EINT for INT------------
#find the intersection points of line segments
ya=scale(INTcomp)   # ya=(ACTcomp+EINTcomp)/sqrt(2) 
yb=ACTcomp/sqrt(2) -  (ACTcomp/sqrt(2))[1] + ya[1] - 0.0001  #0.0001 is added just for plotting technic. 

above=ya>yb
cross=which(diff(above)!=0)
slopes1=ya[cross+1]-ya[cross]
slopes2=yb[cross+1]-yb[cross]
x.cross=cross + ((yb[cross] - ya[cross]) / (slopes1-slopes2))
y.cross=ya[cross] + (slopes1*(x.cross-cross))

# Calc. polygons
x=1:length(year) ; x2=sort(c(x,x.cross))
k=which(x2%%1 != 0)
ya2=rep(NA,length(x2)) ; ya2[k]=y.cross ; ya2[which(is.na(ya2))]=ya
yb2=rep(NA,length(x2)) ; yb2[k]=y.cross ; yb2[which(is.na(yb2))]=yb

#plot positive area
plot(-999, -999, xlim=c(1984, 2015), ylim=c(-3, 4), axes=FALSE, xlab='', ylab='')

k2=which(ya2>=yb2)
k3=c(k2[1],k2[which(diff(k2)!=1)+1])
k4=c(k2[which(diff(k2)!=1)],rev(k2)[1])
for (i in 1:length(k3)) 
{x3=c(x2[k3[i]:k4[i]],rev(x2[k3[i]:k4[i]])) 
 y3=c(ya2[k3[i]:k4[i]],yb2[k4[i]:k3[i]])
  polygon(year[1]-1+x3,y3,border=NA,col=colset[7])
 }  

#plot negative area
k2=which(ya2<=yb2)
k3=c(k2[1],k2[which(diff(k2)!=1)+1])
k4=c(k2[which(diff(k2)!=1)],rev(k2)[1])
for (i in 1:length(k3)) 
{x4=c(x2[k3[i]:k4[i]],rev(x2[k3[i]:k4[i]])) 
  y4=c(ya2[k3[i]:k4[i]],yb2[k4[i]:k3[i]])
  polygon(year[1]-1+x4,y4,border=NA,col=colset[8])
 }  

abline(h=ya[1], lty=3)
lines(year,ya,col=colset[1],lwd=3)
lines(year,yb,col=colset[3],lwd=3)
lines(year,scale(ACEcomp)*sd(yb)+mean(yb),col=3,lwd=2,lty=1)

#axis
axis(1, at = seq(1980, 2020, 5), labels=seq(1980, 2020, 5), cex.axis=1)
intlev=seq(30,120,10);  intsdlev=(intlev-mean(INTjt))/sd(INTjt) 
axis(2, at = intsdlev, labels=intlev, las=1, cex.axis=1, col=1, col.axis=1)
axis(4, at = intsdlev, labels=rep('',length(intsdlev)), las=1, cex.axis=1, col=1, col.axis=1)

mtext('Year', side=1, line=2.6, cex=1.2)
mtext(expression(paste('Mean LMI (m s'^-1,')')), side=2, line=2.5, cex=1.2, col=1)

legend('topleft',c('INT variation','ACT variation starting from 1984 level (scaled to INT)','ACE variation (scaled to ACT)',
'intensification by EINT','weakening by EINT'), # text
bty='n',
pch=c(-1,-1,-1,15,15), # -1 is a nonexistent point
col=c(colset[1],colset[3],3, colset[7],colset[8]), 
pt.cex=1.5,
lty=c(1,1,1,0,0), # sequence of line types for the legend
lwd=c(3,3,2,3,3), cex=0.8)
box()

Fig. 4. WNP TC intensity (INT) compared to the portion of activity (ACT) during JJASON over the 32 years (1984–2015). The gap between INT and ACT is the portion of the efficiency of INT (EINT), which is the intensification at the expense of TC occurrences (blasting effect). ACT is also compared to Accumulated Cyclone Energy (ACE), representing the similarity between the two metrics. TC activity indicated by ACT and ACE are seen to have no trend, while EINT shows an increasing trend.

# Ref.  Probability level of the occurrences of FRQ, INT, ACT and EINT in 2015
##############################################################

Fn=ecdf(FRQcomp);Fn(FRQcomp[32])

## [1] 0.3125

# 0.3125
Fn=ecdf(INTcomp);Fn(INTcomp[32])

## [1] 1

# 1
Fn=ecdf(ACTcomp);Fn(ACTcomp[32])

## [1] 0.71875

# 0.71875
 Fn=ecdf(EINTcomp);Fn(EINTcomp[32])

## [1] 0.96875

# 0.96875

Figure 5. Contribution of warming environment to 2015 EINT

rm(list=ls())

colset=c('black','#ff7300','blue','green','purple','#e1ff00', '#ffb400', '#9933cc30', '#5500ff')

#------------read data-------------------------------------
year=1984:2015

NSOI=-1*read.table('Data/SOIann_jjason_1951to2015.dat',T)[(year[1]-1950):(rev(year)[1]-1950),1]  
GMSST=read.table('Data/GMSSTann_jjason_1854to2015.dat',T)[(year[1]-1853):(rev(year)[1]-1853),1]
sNSOI=scale(NSOI)[,1]
sGMSST=scale(GMSST)[,1]

INTjt=read.table('Data/intjtwc_jjason.txt',T)[,1]*0.51444  #kt to m/s
INTjm=read.table('Data/intjma_jjason.txt',T)[,1]*0.51444  #kt to m/s
FRQjt=read.table('Data/frqjtwc_jjason.txt',T)[,1]
FRQjm=read.table('Data/frqjma_jjason.txt',T)[,1]
 
#composite
INTcomp=((scale(INTjt)+scale(INTjm))/sqrt(2))
FRQcomp=((scale(FRQjt)+scale(FRQjm))/sqrt(2))
EINTcomp=(scale(INTcomp)-scale(FRQcomp))/sqrt(2)

#------------plot  role of GMSST for EINT------------

#ya=scale(INTcomp)   # ya=(ACTcomp+EINTcomp)/sqrt(2) 
#yb=rep(0,32)
#yb=ACTcomp/sqrt(2) -  (ACTcomp/sqrt(2))[1] + ya[1] - 0.0001  #0.0001 is added just for plotting technic. 

EINTp1 = as.numeric(predict(lm(EINTcomp~GMSST)))
EINTp2 = as.numeric(predict(lm(EINTcomp~NSOI +GMSST)))

ya= (EINTcomp-EINTcomp[1])/sqrt(2) 
yb=rep(0,32)+0.0001
yc= (EINTp1-EINTcomp[1])/sqrt(2)
yc2= (EINTp2-EINTcomp[1])/sqrt(2) 

#find the intersection points of line segments
above=ya>yb
cross=which(diff(above)!=0)
slopes1=ya[cross+1]-ya[cross]
slopes2=yb[cross+1]-yb[cross]
x.cross=cross + ((yb[cross] - ya[cross]) / (slopes1-slopes2))
y.cross=ya[cross] + (slopes1*(x.cross-cross))

# Calc. polygons
x=1:length(year) ; x2=sort(c(x,x.cross))
k=which(x2%%1 != 0)
ya2=rep(NA,length(x2)) ; ya2[k]=y.cross ; ya2[which(is.na(ya2))]=ya
yb2=rep(NA,length(x2)) ; yb2[k]=y.cross ; yb2[which(is.na(yb2))]=yb

#plot positive area
plot(-999, -999, xlim=c(1984, 2015), ylim=c(-1, 3), axes=FALSE, xlab='', ylab='')

k2=which(ya2>=yb2)
k3=c(k2[1],k2[which(diff(k2)!=1)+1])
k4=c(k2[which(diff(k2)!=1)],rev(k2)[1])
for (i in 1:length(k3)) 
{x3=c(x2[k3[i]:k4[i]],rev(x2[k3[i]:k4[i]])) 
 y3=c(ya2[k3[i]:k4[i]],yb2[k4[i]:k3[i]])
  polygon(year[1]-1+x3,y3,border=NA,col=colset[7])
 }  

#plot negative area
k2=which(ya2<=yb2)
k3=c(k2[1],k2[which(diff(k2)!=1)+1])
k4=c(k2[which(diff(k2)!=1)],rev(k2)[1])
for (i in 1:length(k3)) 
{x4=c(x2[k3[i]:k4[i]],rev(x2[k3[i]:k4[i]])) 
  y4=c(ya2[k3[i]:k4[i]],yb2[k4[i]:k3[i]])
  polygon(year[1]-1+x4,y4,border=NA,col=colset[8])
 }  

abline(h=0, lty=3)
lines(year,yc, col=colset[2],lwd=4)     # modeled EINT by GMSST
lines(year,yc2,col=colset[9],lwd=2)      # modeled EINT by NSOI + GMSST

#axis
axis(1, at = seq(1980, 2020, 5), labels=seq(1980, 2020, 5), cex.axis=1)
Dintlev=seq(-20,20,5) ;  intsdlev=Dintlev/sd(INTjt)
axis(2, at = intsdlev, labels=Dintlev, las=1, cex.axis=1, col=1, col.axis=1)
axis(4, at = intsdlev, labels=rep('',length(intsdlev)), las=1, cex.axis=1, col=1, col.axis=1)

mtext('Year', side=1, line=2.6, cex=1.2)
mtext(expression(paste(Delta, ' Mean LMI (m s'^-1,')')), side=2, line=2.5, cex=1.2, col=1)

legend('topleft',
    c('intensification by EINT (departure from 1984)', 'weakening by EINT (departure from 1984)',
       'modeled EINT by GMSST',  'modeled EINT by GMSST and NSOI'), # text  
bty='n',
pch=c(15,15, -1,-1), # -1 is a nonexistent point
col=c(colset[7],colset[8],colset[2],colset[9]), 
pt.cex=1.5,
lty=c(0,0,1,1), # sequence of line types for the legend
lwd=c(0,0,4,2), cex=0.8)
box()

Fig. 5. EINT departure from 1984 and its modeled results. All values are averaged for JJASON and represented as INT scale. Regression model of EINT reveals that 2015 EINT is contributed mostly by the warming environment indicated by global mean sea surface temperature (GMSST). Additional contribution of El Nino environment indicated by negative Southern Oscillation Index (NSOI) seems negligible to the highest level of 2015 EINT.

# Ref.  Correlation analysis among variables
##############################################################

S1=cor.test(NSOI,year) ; S1$estimate ; S1$conf.int

##         cor 
## -0.09512011

## [1] -0.4295663  0.2622729
## attr(,"conf.level")
## [1] 0.95

S2=cor.test(NSOI,GMSST)  ; S2$estimate ; S2$conf.int

##      cor 
## 0.161862

## [1] -0.1980076  0.4832795
## attr(,"conf.level")
## [1] 0.95

S3=cor.test(EINTcomp,predict(lm(EINTcomp~GMSST)))  ; S3$estimate ; S3$conf.int 

##       cor 
## 0.7124762

## [1] 0.4840336 0.8499995
## attr(,"conf.level")
## [1] 0.95

S4=cor.test(EINTcomp,predict(lm(EINTcomp~GMSST+NSOI)))  ; S4$estimate ; S4$conf.int

##       cor 
## 0.7267025

## [1] 0.5063045 0.8579857
## attr(,"conf.level")
## [1] 0.95

100*(round((S3$estimate^2),3))

##  cor 
## 50.8

100*(round((S4$estimate^2 -S3$estimate^2),3))

## cor 
##   2

[Supplementary] Profiles of MSE and GPH in 2015

rm(list=ls())

par(mfrow=c(1,1),mar=c(5,6,3,4))

#------------read data-------------------------------------
year=1984:2015

NSOI=-1*read.table('Data/SOIann_jjason_1951to2015.dat',T)[(year[1]-1950):(rev(year)[1]-1950),1]  
GMSST=read.table('Data/gmsst_jjason_1854to2015.dat',T)[(year[1]-1853):(rev(year)[1]-1853),1]

z=read.table('Data/hgt.1984to2015.jjason.wnp.txt',FALSE)
q=read.table('Data/shum.1984to2015.jjason.wnp.txt',FALSE)
t=read.table('Data/temp.1984to2015.jjason.wnp.txt',FALSE)

#-------------calculate annual MSE------------------------------------
Cp= 1004 # specific heat capacity of air (J/(kg.K))
L= 2.5*10^6 # latent heat of vaporization (J/kg)
g= 9.81 # accerelation due to gravity (m/s2)

SHann=Cp*(t+273.15)        # sensible heating energy (J/kg)
LHann=L*q/1000                  # latent heating (J/kg)  # q over 300hPa is ignored to 0
PEann=g*z                            # potential energy (m2/s2)=(J/kg)
MSEann=SHann+PEann+LHann              #annual Moist static energy

#--------make each profile (30S~30N area-weighted mean values)---
Year=1984:2015
Lat=seq(-90,90,2.5)
n1=which(Lat==0) ;  n2=which(Lat==30)   # tropics

var=t
profileset=NULL       #dim[time,lev]
 for (lev in 1:17)              # for each lev (1000 hPa to 10 hPa)
 {tmsrs=NULL
  annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
  for (ann in annualnumset)   # time-series of area weight mean (0~30N)
  {area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
   tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}  
   profileset=rbind(profileset,tmsrs)
 }
profileset_t=profileset

var=q
profileset=NULL             #dim[time,lev]
 for (lev in 1:8)              # for each lev (1000 hPa to 10 hPa)
 {tmsrs=NULL
  annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
  for (ann in annualnumset)   # time-series of area weight mean (0~30N)
  {area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
   tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}  
   profileset=rbind(profileset,tmsrs)
 }
profileset_q=profileset

var=z
profileset=NULL             #dim[time,lev]
 for (lev in 1:17)              # for each lev (1000 hPa to 10 hPa)
 {tmsrs=NULL
  annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
  for (ann in annualnumset)   # time-series of area weight mean (0~30N)
  {area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
   tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}  
   profileset=rbind(profileset,tmsrs)
 }
profileset_z=profileset

var=MSEann
profileset=NULL             #dim[time,lev]
 for (lev in 1:17)              # for each lev (1000 hPa to 10 hPa)
 {tmsrs=NULL
  annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
  for (ann in annualnumset)   # time-series of area weight mean (0~30N)
  {area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
   tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}  
   profileset=rbind(profileset,tmsrs)
 }
profileset_MSE=profileset

#------------s.d.-------------------
sd_MSE=c()
sd_z=c()
for (i in 1:17)
{sd_MSE=rbind(sd_MSE,scale(profileset_MSE[i,])[,1])
  sd_z=rbind(sd_z,scale(profileset_z[i,])[,1])
}

#---------------------------plot-----------------------------------------------
p=c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10)
p2=log(c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10))

plot(-999,-999,xlim=c(-4,4),ylim=c(log(1000),log(150)),xlab='',ylab='',axes=FALSE)
abline(v=0,lty=2,col='gray')

lines(sd_z[,32],p2,col='blue',lwd=5)                #s.d. of z (2015)
lines(sd_MSE[,32],p2,col='orange',lwd=5)              #s.d. of MSE (2015)

axis(1,at=seq(-4,4,1),labels=seq(-4,4,1),cex.axis=1)
axis(2,at=p2,labels=p,cex.axis=1,las=1)  
axis(3,at =seq(-4,4,1),labels=rep('',length=length(seq(-4,4,1))))
axis(4, at=p2,labels=rep('',length(p)))
mtext('Standardized value (s.d.)',side=1,line=3, cex=1.4)
mtext('Pressure (hPa)',side=2,line=3.5, cex=1.4)
legend(-4,log(200), c('MSE 2015','Geopotential height 2015')
,bty='n',col=c('orange','blue'),lwd=2,cex=1)
box()

[Supplementary] Profiles of standardized moist static energy (MSE) and geopotential height for the year 2015 in the tropical region (0–30N) of the western North Pacific (100–180E). All values are averaged for JJASON over 32 years (1984–2015), and then standardized values of 2015 are represented.

END