Global tropical cyclone climate has been investigated with indicators of frequency, intensity and activity. However, a full understanding of global warming's influence on tropical cyclone climate remains elusive because of the incomplete nature of these indicators. Here we form a complete three-dimensional variability space of tropical cyclone climate where the variabilities are continuously linked and find that global ocean warmth best explains the out-of-phase relationship between intensity and frequency of global tropical cyclones. In a year with greater ocean warmth, the tropical troposphere is capped by higher pressure anomaly in the middle and upper troposphere even with higher moist static energy anomaly in the lower troposphere, which is thought to inhibit overall tropical cyclone occurrences but lead to greater intensities. Statistical consequence is represented as the trade-off between intensity and frequency. We calculate an average increase in global tropical cyclone intensity of 1.3 ms\( ^{-1} \) over the past 30 years of ocean warming occurring at the expense of 6.1 tropical cyclones worldwide.
Tropical cyclones (TCs) are perhaps the least welcomed natural phenomena on our planet. Even well developed and highly complex societies are vulnerable to their destructiveness because of significant exposure. Ocean warmth is increasing with global warming and a major concern is how TC climate will change. Recent consensus from theoretical and numerical projections of future TC climate reveals decreasing frequency and increasing intensity in this warming environment. Nevertheless, there is little observational support except for the increasing intensity of the strongest portion of TCs. Trend analysis is a popular approach, but results are constrained by a focus on just a few TC climate indicators (e.g., Accumulated Cyclone Energy). Is it possible that the influence of global warming on TCs is aligned somewhere in-between these indicators?
Here we examine a globally consistent TC response to ocean warming by linearly linking frequency, intensity, and activity in a continuous variability space\( ^{10} \). The linear approach is expected to provide a convenient framework for understanding the TC climate structure. TCs whose lifetime-maximum wind (LMW) speed exceeds 17 ms\( ^{-1} \) are selected annually. In this framework, the annual number of TC occurrences and the annual mean LMW are used as indicators for frequency (F) and intensity (I), respectively. These two primary TC indicators are not orthogonal to each other, but the corresponding eigenvectors derived from a principal component analysis create a continuous TC variability space (Supplementary Fig. 1). The principal component of the in-phase relationship between F and I indicates a variability direction close to those of Accumulated Cyclone Energy (ACE) and Power Dissipation Index (PDI) (Supplementary Fig. 6). This variability direction is named A to stand for 'activity' and note the other principal component indicates a variability direction orthogonal to A. This variability represents the out-of-phase relationship between F and I and is named E to stand for 'efficiency of intensity', since it shows how much more I contributes to A than does F. Practically, E reveals the degree of trade-off between I and F at a given A. If it were not for E, the TC climate variability space with I, F and A would not span the two-dimensional variability space.
The same orthogonal structure is applied to our environmental framework, where the Southern Oscillation Index (SOI) and global sea surface temperature (SST) are chosen as the primary indicators (Supplementary Fig. 3). The former indicates El Nino (N) and the latter global ocean warmth (O). This environmental space suggests two orthogonal variabilities of warm-year El Nino (W) and cold-year El Nino (C), which are indicated by the two principal components. Here, -W and -C point to cold-year La Nina and warm-year La Nina inversely. Note that 'warm' and 'cold' refer to global temperature rather than the El Nino phenomenon itself. Figure 1 demonstrates how the two orthogonal spaces, one defined by TC climate and the other by the TC environment, are placed in a three-dimension variability space. Projection of the TC climate framework onto the environmental framework points to the best explanatory environmental variability for each TC climate variability direction (Supplementary Fig. 4). A projection length (correlation) shows how much the TC climate variability is explained. The contribution of El Nino and global ocean warmth to the projection length is also identified.
par(mar = c(0, 0, 0, 0))
require(png)
## Loading required package: png
require(grid)
## Loading required package: grid
img <- readPNG("figure/Main_Figs/Main_Fig1.png")
grid.raster(img)
Fig. 1 Schematic of a three-dimensional variability space. The TC climate variability space and the environmental variability space are constructed separately. The two variability spaces have a crossing angle and different magnitudes. The projection of TC climate onto the environmental climate shows which environmental variable best explains each TC climate variability direction.
Figure 2 presents the projection results for the globe and six regional TC climates. The past 29-year (1984~2012) observations are used here on annual basis. Each indicator is the proxy of the synthetic environmental structure and is circularly chained. A and E have fixed positions since they are orthogonal to each other by definition, whereas F and I position varies depending on TC climate (Supplementary Fig. 2). The best explanatory environment for each TC climate is marked in the upper abscissa. The response of western North Pacific TCs to the defined environment (both El Ni~{n}o and global warming) is noticeably high and significant (\( \alpha\leq 0.05 \)) at all directions (red line). Significant response in the eastern North Pacific and the North Atlantic occurs in a selective variability range. TCs in the North Indian and the Southern Hemisphere are not well explained by this framework. This is not due to lower statistical power since the number of TCs is not the sample size. Each basin has 29 annual values. Nonetheless, the global TC climate coming from these regional features is seen to have a significant relationship with the defined environment at all directions.
It is found that El Nino is the best explanatory variability for global TC activity. On the other hand, global ocean warmth is the best explanatory variability for the efficiency of intensity as indicated by the vertical green lines in Fig. 2. By construction, the efficiency of intensity is the intensity variation with activity variation removed. This result implies that global ocean warmth increases global TC intensity but decreases frequency at the same time. Hence, the efficiency level explained by global ocean warmth is the amount (magnitude) of trade-off between intensity and frequency. Regionally the western North Pacific TC climate variability is similar to that of the globe as a whole. There is no direct link between the two since regional values of intensity, activity, and efficiency are determined by the regional frequency.
rm(list = ls())
require(akima) # for interp --> interpolation
## Loading required package: akima
require(RColorBrewer) # for Brewer.pal --> color
## Loading required package: RColorBrewer
1.set option
Year = 1984:2012
2.read data
frqgl = read.table("Data/frqgl.txt", TRUE)[, 1]
intgl = read.table("Data/intgl.txt", T)[, 1]
frqni = read.table("Data/frqni.txt", TRUE)[, 1]
intni = read.table("Data/intni.txt", T)[, 1]
frqwnp = read.table("Data/frqwnp.txt", TRUE)[, 1]
intwnp = read.table("Data/intwnp.txt", T)[, 1]
frqenp = read.table("Data/frqenp.txt", TRUE)[, 1]
intenp = read.table("Data/intenp.txt", T)[, 1]
frqna = read.table("Data/frqna.txt", TRUE)[, 1]
intna = read.table("Data/intna.txt", T)[, 1]
frqsi = read.table("Data/frqsi.txt", TRUE)[, 1]
intsi = read.table("Data/intsi.txt", T)[, 1]
frqsp = read.table("Data/frqsp.txt", TRUE)[, 1]
intsp = read.table("Data/intsp.txt", T)[, 1]
gmsst = read.table("Data/gmsst_12m_1854to2013.dat", TRUE)[(Year[1] - 1853):(rev(Year)[1] -
1853), 1]
nsoi = -1 * read.table("Data/SOIann_12m_1951to2013.dat", TRUE)[(Year[1] - 1950):(rev(Year)[1] -
1950), 1]
3.global result
par(mfrow = c(1, 1), mar = c(5, 6, 6, 6), bg = "#333333", col.axis = "white",
col.lab = "white")
frq = frqgl
int = intgl
# extract orthogonal vectors
act = (scale(int) + scale(frq))/sqrt(2) # Activity
eint = (scale(int) - scale(frq))/sqrt(2) # Efficiency of INTensity
# find the best explanatory Env indicator for each TC climate indicator
# --angles between primary indicators in their orthogonal frameworks
Z = acos(cor(frq, int)) * 180/pi # angle between frq and int
A = acos(cor(gmsst, nsoi)) * 180/pi # angle between gmsst & nsoi
# --variability chain
TCloc = c(-90, -Z/2, 0, Z/2, 90, 180 - Z/2, 180, 180 + Z/2)
# -E,F,A, I, E, -F,-A,-I (when setting A to 0)
ENVloc = c(-90, -A/2, 0, A/2, 90, 180 - A/2, 180, 180 + A/2)
# -C,O, W,N, C, -O, -W,-N (when setting W to 0)
# --find the best explanatory Env variability directions
TCang = seq(-90, (270 - 5), 5) # TC angle--> TC indicator
ENVang = NULL # ENV angle--> Env indicator
for (deg in TCang) {
TCindicator = cos(deg * pi/180) * scale(act) + sin(deg * pi/180) * scale(eint)
best = predict(lm(TCindicator ~ gmsst + nsoi))
cos_best = cor(TCindicator, best) # refer to Suppl.Fig.4
cos_nsoi = cor(TCindicator, nsoi) # refer to Suppl.Fig.4
cos_gmsst = cor(TCindicator, gmsst) # refer to Suppl.Fig.4
B = acos(cos_nsoi/cos_best) * 180/pi
C = acos(cos_gmsst/cos_best) * 180/pi
if (B <= A & C <= A) {
D = A/2 - B
if (D < 0)
D = 360 + D
}
if (B > A | C > A) {
F = 1
if (C > B)
D = B + A/2
if (C < B)
D = 360 - (C + A/2)
}
if (B > (180 - A) & C > (180 - A))
D = B + A/2
ENVang = c(ENVang, D)
}
# get the TC climate variabilities matching with the Env chain /This is for
# plotting the axis marks /predict from spline result
x = ENVang
y = TCang
x2 = x[order(x)]
y2 = y[order(x)]
k = which.max(abs(diff(y2)))
if (y2[k] == -90) {
y2[1:k] = y2[1:k] + 360
} else {
y2[1:k] = y2[1:k] - 360
}
x3 = c(x2 - 360, x2, x2 + 360)
if (diff(y2)[1] < 0) y3 = c(y2 + 360, y2, y2 - 360)
if (diff(y2)[1] > 0) y3 = c(y2 - 360, y2, y2 + 360)
interp = smooth.spline(x3, y3, spar = 0)
pred = predict(interp, ENVloc[1:4])
xdeg = pred$x
ydeg = pred$y
imsi = rep(NA, 8)
imsi[1:4] = ydeg[1:4]
imsi[5:8] = ydeg[1:4] + 180
imsi = imsi%%360
imsi[imsi > 270] = imsi[imsi > 270] - 360
ENVonTCang = round(imsi, 1)
# compute the correlations and p-values
TCang = seq(-90, 270, 2.5) # TC angle
corr_TandBEST = NULL
pval_total = NULL
corr_BESTandO = NULL
corr_BESTandN = NULL
for (deg in TCang) {
T = cos(deg * pi/180) * scale(act) + sin(deg * pi/180) * scale(eint)
BEST = predict(lm(T ~ gmsst + nsoi))
ctest = cor.test(T, BEST)
corr_TandBEST = c(corr_TandBEST, ctest$estimate)
pval_total = c(pval_total, ctest$p.val)
corr_BESTandO = c(corr_BESTandO, cor(BEST, gmsst)) # O: gmsst
corr_BESTandN = c(corr_BESTandN, cor(BEST, nsoi)) # N: nsoi
}
r.sq = 100 * (corr_TandBEST^2) # proportion that BEST can explain T
r.sq_O = r.sq * (corr_BESTandO^2) # O portion
r.sq_N = r.sq * (corr_BESTandN^2) # N portion
# plot
# This figure shows how much env indicator explains TC climate. Besides,
# the contributions of global ocean warmth and El Nino are shown.
plot(-999, -999, xlim = c(-90, 270), ylim = c(0, 50), xlab = "", ylab = "",
axes = FALSE, xaxs = "i", yaxs = "i")
polygon(c(-180, 360, 360, -180), c(-50, -50, 150, 150), border = NA, col = "white")
# --for portion_gmsst
xlist = c(TCang - 360, TCang, TCang + 360)
ylist1 = rep(r.sq, 3)
ylist2 = rep(r.sq - r.sq_O, 3)
polygon(c(xlist, rev(xlist)), c(ylist1, rev(ylist2)), border = NA, col = "#ffff0099")
# --for portion_nsoi
xlist = c(TCang - 360, TCang, TCang + 360) # repeat to fill each sides of the figure box
ylist1 = rep(r.sq_N, 3)
polygon(c(xlist, rev(xlist)), c(ylist1, rep(0, length(xlist))), border = NA,
col = "#99339988")
# --for total variance
imsi = rep(NA, length(pval_total))
imsi[pval_total <= 0.05] = r.sq[pval_total <= 0.05]
ylist = rep(imsi, 3)
xlist = c(TCang - 360, TCang, TCang + 360)
lines(xlist, ylist, lwd = 2, col = "#ff0000")
# --line for GMSST
abline(v = ENVonTCang[c(2, 6)], col = "#66ff00", lwd = 2)
# --axis
axis(1, at = c(TCloc, 270), labels = c("-E", "F", "A", "I", "E", "-F", "-A",
"-I", "-E"), cex.axis = 1.5, col = "white")
axis(2, at = seq(0, 50, 10), labels = seq(0, 50, 10), cex.axis = 1.5, las = 1,
col = "white")
axis(3, at = ENVonTCang, labels = c("-C", "O", "W", "N", "C", "-O", "-W", "-N"),
cex.axis = 1.5, col = "white")
axis(4, at = 100 * (seq(0, 1, 0.1)^2), labels = c(0, "", seq(0.2, 1, 0.1)),
cex.axis = 1.5, las = 1, col = "white")
mtext("TC Climate", side = 1, line = 3.8, cex = 2.2, col = "white")
mtext("Best Explanatory Environment", side = 3, line = 3.8, cex = 2.2, col = "white")
mtext("Variance Portion (%)", side = 2, line = 3.5, cex = 2, col = "white")
text(par("usr")[2] + 44, 42, srt = -90, adj = 0, labels = "Correlation Coefficient",
cex = 1.8, xpd = TRUE, col = "white")
text(-94, 47, adj = 0, " Global", cex = 2.5, font = 1)
mtext("a", side = 3, adj = -0.1, line = 3.6, cex = 3.2, col = "white", font = 2)
4.note on variability chains
par(mfrow = c(1, 1), mar = c(5, 6, 6, 6), bg = "#333333", col.axis = "white",
col.lab = "white")
plot(-999, -999, xlim = c(0, 100), ylim = c(0, 100), xlab = "", ylab = "", axes = FALSE)
text(0, 100, "Environmental Variability Chain", col = "white", cex = 1.3, font = 1.2,
adj = 0)
segments(0, 93, 90, 93, col = "white", lwd = 2)
hh = seq(86, 60, length = 4)
text(rep(0, 4), hh, c("C", "N", "W", "O"), col = "white", cex = 1.5, font = 2,
adj = 0)
Name = bquote(cold - year * " " * El * " " * Ni * tilde(n) * o)
text(10, hh[1], Name, col = "white", cex = 1.3, adj = 0)
Name = bquote(El * " " * Ni * tilde(n) * o)
text(10, hh[2], Name, col = "white", cex = 1.3, adj = 0)
Name = bquote(warm - year * " " * El * " " * Ni * tilde(n) * o)
text(10, hh[3], Name, col = "white", cex = 1.3, adj = 0)
text(10, hh[4], "global ocean warmth", col = "white", cex = 1.3, adj = 0)
k = 60
text(0, 100 - k, "TC Climate Variability Chain", col = "white", cex = 1.3, font = 1,
adj = 0)
segments(0, 93 - k, 90, 93 - k, col = "white", lwd = 2)
hh = seq(86 - k, 60 - k, length = 4)
text(rep(0, 4), hh, c("F", "A", "I", "E"), col = "white", cex = 1.5, font = 2,
adj = 0)
text(rep(10, 4), hh[1:4], c("frequency", "activity", "intensity", "efficiency of intensity"),
col = "white", cex = 1.3, adj = 0)
5.loop for regions
par(mfrow = c(1, 1), mar = c(5, 6, 6, 6), bg = "#333333", col.axis = "white",
col.lab = "white")
for (bsn in 2:7) {
if (bsn == 2) {
frq = frqwnp
int = intwnp
}
if (bsn == 3) {
frq = frqna
int = intna
}
if (bsn == 4) {
frq = frqenp
int = intenp
}
if (bsn == 5) {
frq = frqni
int = intni
}
if (bsn == 6) {
frq = frqsi
int = intsi
}
if (bsn == 7) {
frq = frqsp
int = intsp
}
# extract orthogonal vectors
act = (scale(int) + scale(frq))/sqrt(2) # Activity
eint = (scale(int) - scale(frq))/sqrt(2) # Efficiency of INTensity
# find the best explanatory Env indicator for each TC climate indicator
# --angles between primary indicators in their orthogonal frameworks
Z = acos(cor(frq, int)) * 180/pi # angle between frq and int
A = acos(cor(gmsst, nsoi)) * 180/pi # angle between gmsst & nsoi
# --variability chain
TCloc = c(-90, -Z/2, 0, Z/2, 90, 180 - Z/2, 180, 180 + Z/2)
# -E,F,A, I, E, -F,-A,-I (when setting A to 0)
ENVloc = c(-90, -A/2, 0, A/2, 90, 180 - A/2, 180, 180 + A/2)
# -C,O, W,N, C, -O, -W,-N (when setting W to 0)
# --find the best explanatory Env variability directions
TCang = seq(-90, (270 - 5), 5) # TC angle--> TC indicator
ENVang = NULL # ENV angle--> Env indicator
for (deg in TCang) {
TCindicator = cos(deg * pi/180) * scale(act) + sin(deg * pi/180) * scale(eint)
best = predict(lm(TCindicator ~ gmsst + nsoi))
cos_best = cor(TCindicator, best) # refer to Suppl.Fig.4
cos_nsoi = cor(TCindicator, nsoi) # refer to Suppl.Fig.4
cos_gmsst = cor(TCindicator, gmsst) # refer to Suppl.Fig.4
B = acos(cos_nsoi/cos_best) * 180/pi
C = acos(cos_gmsst/cos_best) * 180/pi
if (B <= A & C <= A) {
D = A/2 - B
if (D < 0)
D = 360 + D
}
if (B > A | C > A) {
F = 1
if (C > B)
D = B + A/2
if (C < B)
D = 360 - (C + A/2)
}
if (B > (180 - A) & C > (180 - A))
D = B + A/2
ENVang = c(ENVang, D)
}
# get the TC climate variabilities matching with the Env chain /This is for
# plotting the axis marks /predict from spline result
x = ENVang
y = TCang
x2 = x[order(x)]
y2 = y[order(x)]
k = which.max(abs(diff(y2)))
if (y2[k] == -90) {
y2[1:k] = y2[1:k] + 360
} else {
y2[1:k] = y2[1:k] - 360
}
x3 = c(x2 - 360, x2, x2 + 360)
if (diff(y2)[1] < 0)
y3 = c(y2 + 360, y2, y2 - 360)
if (diff(y2)[1] > 0)
y3 = c(y2 - 360, y2, y2 + 360)
interp = smooth.spline(x3, y3, spar = 0)
pred = predict(interp, ENVloc[1:4])
xdeg = pred$x
ydeg = pred$y
imsi = rep(NA, 8)
imsi[1:4] = ydeg[1:4]
imsi[5:8] = ydeg[1:4] + 180
imsi = imsi%%360
imsi[imsi > 270] = imsi[imsi > 270] - 360
ENVonTCang = round(imsi, 1)
# compute the correlations and p-values
TCang = seq(-90, 270, 2.5) # TC angle
corr_TandBEST = NULL
pval_total = NULL
corr_BESTandO = NULL
corr_BESTandN = NULL
for (deg in TCang) {
T = cos(deg * pi/180) * scale(act) + sin(deg * pi/180) * scale(eint)
BEST = predict(lm(T ~ gmsst + nsoi))
ctest = cor.test(T, BEST)
corr_TandBEST = c(corr_TandBEST, ctest$estimate)
pval_total = c(pval_total, ctest$p.val)
corr_BESTandO = c(corr_BESTandO, cor(BEST, gmsst)) # O: gmsst
corr_BESTandN = c(corr_BESTandN, cor(BEST, nsoi)) # N: nsoi
}
r.sq = 100 * (corr_TandBEST^2) # proportion that BEST can explain T
r.sq_O = r.sq * (corr_BESTandO^2) # O portion
r.sq_N = r.sq * (corr_BESTandN^2) # N portion
# plot
# This figure shows how much env indicator explains TC climate. Besides,
# the portions of GMSST and NSOI are shown.
plot(-999, -999, xlim = c(-90, 270), ylim = c(0, 50), xlab = "", ylab = "",
axes = FALSE, xaxs = "i", yaxs = "i")
polygon(c(-180, 360, 360, -180), c(-50, -50, 150, 150), border = NA, col = "white")
# --for portion_gmsst
xlist = c(TCang - 360, TCang, TCang + 360)
ylist1 = rep(r.sq, 3)
ylist2 = rep(r.sq - r.sq_O, 3)
polygon(c(xlist, rev(xlist)), c(ylist1, rev(ylist2)), border = NA, col = "#ffff0099")
# --for portion_nsoi
xlist = c(TCang - 360, TCang, TCang + 360) # repeat to fill each sides of the figure box
ylist1 = rep(r.sq_N, 3)
polygon(c(xlist, rev(xlist)), c(ylist1, rep(0, length(xlist))), border = NA,
col = "#99339988")
# --for total variance
imsi = rep(NA, length(pval_total))
imsi[pval_total <= 0.05] = r.sq[pval_total <= 0.05]
ylist = rep(imsi, 3)
xlist = c(TCang - 360, TCang, TCang + 360)
lines(xlist, ylist, lwd = 2, col = "#ff0000")
# --line for GMSST
abline(v = ENVonTCang[c(2, 6)], col = "#66ff00", lwd = 2)
# --axis
axis(1, at = c(TCloc, 270), labels = c("-E", "F", "A", "I", "E", "-F", "-A",
"-I", "-E"), cex.axis = 1.5, col = "white")
axis(2, at = seq(0, 50, 10), labels = seq(0, 50, 10), cex.axis = 1.5, las = 1,
col = "white")
axis(3, at = ENVonTCang, labels = c("-C", "O", "W", "N", "C", "-O", "-W",
"-N"), cex.axis = 1.5, col = "white")
axis(4, at = 100 * (seq(0, 1, 0.1)^2), labels = c(0, "", seq(0.2, 1, 0.1)),
cex.axis = 1.5, las = 1, col = "white")
mtext("TC Climate", side = 1, line = 3.8, cex = 2.2, col = "white")
mtext("Best Explanatory Environment", side = 3, line = 3.8, cex = 2.2, col = "white")
mtext("Variance Portion (%)", side = 2, line = 3.5, cex = 2, col = "white")
text(par("usr")[2] + 44, 42, srt = -90, adj = 0, labels = "Correlation Coefficient",
cex = 1.8, xpd = TRUE, col = "white")
if (bsn == 2)
text(-94, 47, adj = 0, " western North Pacific", cex = 2, font = 1)
if (bsn == 2)
mtext("b1", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
if (bsn == 3)
text(-94, 47, adj = 0, " North Atlantic", cex = 2, font = 1)
if (bsn == 3)
mtext("b2", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
if (bsn == 4)
text(-94, 47, adj = 0, " eastern North Pacific", cex = 2, font = 1)
if (bsn == 4)
mtext("b3", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
if (bsn == 5)
text(-94, 47, adj = 0, " North Indian", cex = 2, font = 1)
if (bsn == 5)
mtext("b4", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
if (bsn == 6)
text(-94, 47, adj = 0, " South Indian", cex = 2, font = 1)
if (bsn == 6)
mtext("b5", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
if (bsn == 7)
text(-94, 47, adj = 0, " South Pacific", cex = 2, font = 1)
if (bsn == 7)
mtext("b6", side = 3, adj = -0.1, line = 3.4, cex = 3.2, col = "white",
font = 2)
}
Fig. 2 Projection length of (a) global and (b) regional TC climate variabilities. Variance portion (%) shows how much the best explanatory environment (upper abscissa) explains each TC climate variability (lower abscissa). Significance of the correlation coefficient is examined using t-values. The ranges of significant correlation (\( \alpha\leq 0.05 \)) are shown with red lines. The contribution of El Nino (N) is shaded in purple and the contribution from global ocean warmth (O) is shaded in yellow. The two have overlapping area which means the portion of positive collinearity. No color is given to the portion of negative collinearity. Green vertical lines indicate the variability direction of O and -O.
What physical mechanisms create the large correlation between efficiency of intensification and ocean warming (large yellow area above E)? We pay special attention to the moisture increases in the tropics (30S-30N). By the Clausius-Clapeyron relationship, more moisture is evaporated from a warmer ocean surface resulting in water vapor concentrating in the lower troposphere. Figure 3a is a correlation map that shows the vertical structure of moist static energy change by different environmental variability directions. Compared with other large-scale environmental changes, global ocean warming is clearly seen to be accompanied by a convectively more unstable troposphere. Moist air pumped from the lower troposphere moves upward along the moist adiabat by which the upper air is heated.
In contrast to this warm anomaly prevailing in the middle and the upper troposphere in the tropics, moisture itself in the lower troposphere is not effectively transported out to the free atmosphere (Supplementary Fig. 10). The consequences are warmer and drier conditions aloft, which is well reflected in the high pressure anomaly over the middle and the upper troposphere (Fig. 3b). With the smaller Coriolis effect at lower latitudes, the high pressure anomaly is well stratified and dominant over the tropics. Thus the two collinear environments, which are induced by the lower tropospheric moisture increase in association with the warmer global ocean, can be characterized as a more unstable troposphere but stronger high pressure in the middle and the upper troposphere. The high pressure anomaly concentrates convection in time and space. These thermodynamic conditions over the global tropical cyclone activity imply fewer TCs, but efficient intensification once a TC is spawned to maintain a constant activity level. Thus TC intensity increases at the expense of TC frequency. Additionally, a dynamic factor is examined using Hadley circulation strength (HCS). HCS is defined by average of the maximum absolute mass streamfunction in each hemisphere. Its connection with the global ocean warmth, however, appears insignificant at least over the period studied here (Supplementary. Fig. 11). This study is expected to provide a good basis for further understanding the mass flux response to global warming, which has been mostly examined by numerical models.
require(akima) # for interp --> interpolation
require(fields) # for image.plot
## Loading required package: fields
## Loading required package: spam
## Spam version 0.41-0 (2014-02-26) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
##
## Loading required package: maps
require(RColorBrewer) # for Brewer.pal --> color
1.set option
Year = 1984:2012
2.read data
frq = read.table("Data/frqgl.txt", TRUE)[, 1]
int = read.table("Data/intgl.txt", TRUE)[, 1]
t = read.table("Data/temp.1984to2013.txt", FALSE) # temperature (K)-273.15
q = read.table("Data/shum.1984to2013.txt", FALSE) # specific humidity (g/kg)
z = read.table("Data/hgt.1984to2013.txt", FALSE) # geopotential height (m)
gmsst = read.table("Data/gmsst_12m_1854to2013.dat", TRUE)[(Year[1] - 1853):(rev(Year)[1] -
1853), 1]
nsoi = -1 * read.table("Data/SOIann_12m_1951to2013.dat", TRUE)[(Year[1] - 1950):(rev(Year)[1] -
1950), 1]
3.find the angles for plotting matching axis marks
# angles between primary indicators in their orthogonal frameworks
Z = acos(cor(frq, int)) * 180/pi # angle between frq and int
A = acos(cor(gmsst, nsoi)) * 180/pi # angle between gmsst & nsoi
# variability chain
TCloc = c(-90, -Z/2, 0, Z/2, 90, 180 - Z/2, 180, 180 + Z/2)
# -E,F,A, I, E, -F,-A,-I (when setting A to 0)
ENVloc = c(-90, -A/2, 0, A/2, 90, 180 - A/2, 180, 180 + A/2)
# -C,O, W,N, C, -O, -W,-N (when setting W to 0)
# find the best explanatory Env variability directions
act = (scale(int) + scale(frq))/sqrt(2) # Activity
eint = (scale(int) - scale(frq))/sqrt(2) # Efficiency of INTensity
TCang = seq(-90, (270 - 5), 5) # TC angle--> TC indicator
ENVang = NULL # ENV angle--> Env indicator
for (deg in TCang) {
TCindicator = cos(deg * pi/180) * scale(act) + sin(deg * pi/180) * scale(eint)
best = predict(lm(TCindicator ~ gmsst + nsoi))
cos_best = cor(TCindicator, best) # refer to Suppl.Fig.4
cos_nsoi = cor(TCindicator, nsoi) # refer to Suppl.Fig.4
cos_gmsst = cor(TCindicator, gmsst) # refer to Suppl.Fig.4
B = acos(cos_nsoi/cos_best) * 180/pi
C = acos(cos_gmsst/cos_best) * 180/pi
if (B <= A & C <= A) {
D = A/2 - B
if (D < 0)
D = 360 + D
}
if (B > A | C > A) {
F = 1
if (C > B)
D = B + A/2
if (C < B)
D = 360 - (C + A/2)
}
if (B > (180 - A) & C > (180 - A))
D = B + A/2
ENVang = c(ENVang, D)
}
# get the TC climate variabilities matching with the Env chain /This is for
# plotting the axis marks /predict from spline result
x = ENVang
y = TCang
x2 = x[order(x)]
y2 = y[order(x)]
k = which.max(abs(diff(y2)))
if (y2[k] == -90) {
y2[1:k] = y2[1:k] + 360
} else {
y2[1:k] = y2[1:k] - 360
}
x3 = c(x2 - 360, x2, x2 + 360)
if (diff(y2)[1] < 0) y3 = c(y2 + 360, y2, y2 - 360)
if (diff(y2)[1] > 0) y3 = c(y2 - 360, y2, y2 + 360)
interp = smooth.spline(x3, y3, spar = 0)
pred = predict(interp, ENVloc[1:4])
xdeg = pred$x
ydeg = pred$y
imsi = rep(NA, 8)
imsi[1:4] = ydeg[1:4]
imsi[5:8] = ydeg[1:4] + 180
imsi = imsi%%360
imsi[imsi > 270] = imsi[imsi > 270] - 360
ENVonTCang = round(imsi, 1)
4.calculate annual MSE (Moist Static Energy)
# constants
Cp = 1004 # specific heat capacity of air (J/(kg.K))
L = LHV = 2.5 * 10^6 # latent heat of vaporization (J/kg)
g = 9.81 # accerelation due to gravity (m/s2)
# MSE annual
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
5.make each profile set and plot
par(mar=c(5,8,8,2.6))
for (varnum in 1:2)
{
# set variable and levels
if (varnum==1) var=MSEann
if (varnum==2) var=z
levend= 11 #150 hPa
#c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10)
Lat=seq(-90,90,2.5)
n1=which(Lat==-30) ; n2=which(Lat==30) # tropics
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 (30S~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)
}
# calculate correlations between var and each best explanatory env
degree=seq(-90,270,5) # TC indicator direction (starting from -EINT)
rad=degree*pi/180
corrset = matrix(data = 0, nrow = 17, ncol = length(rad)) # correlation
pvset = matrix(data = 0, nrow = 17, ncol = length(rad)) # two-sided p-value
for (d in 1:length(rad)) # for each indicator
{TCindicator=(cos(rad[d])*scale(act)+sin(rad[d])*scale(eint))/sqrt(2)
pred=predict(lm(TCindicator~gmsst+nsoi)) # this indicates best env variability
for (lev in 1:levend)
{test=cor.test(pred,profileset[lev,])
corrset[lev,d]=test$estimate
pvset[lev,d]=test$p.val
}
}
# plot
#p=c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10)
p=c(1000,925,850,700,600,500,400,300,250,200,150)
Lev=p
#Lev=log(p) # when using log scale p
collist_negative=colorRampPalette(brewer.pal(9,'PuBu'))(10) # extract 10 from Blue-Purple
collist_positive=colorRampPalette(brewer.pal(9,'YlOrRd'))(10) # extract 10 from Yellow-Orange-Red
collist=c(rev(collist_negative),collist_positive) # comprise 20 color set
A1=rep(degree,length(Lev)) # Lat column
A2=rep(Lev,each=length(degree)) # Lev column
A3=NULL; for (j in 1:length(Lev)) A3=c(A3,corrset[j,]) # two sided p-value
INTERP=interp(A1,A2,A3,xo=seq(min(A1),max(A1),length=350),yo=seq(min(A2),max(A2),length=350),duplicate='mean')
image.plot(INTERP,axes=FALSE,col=collist,legend.args=list( text='r',side=3,cex=1.3),xlim=c(-90,270),ylim=c(1000,150),nlevel=20)
A4=NULL; for (j in 1:length(Lev)) A4=c(A4,pvset[j,]) # correlation
INTERP2=interp(A1,A2,A4,xo=seq(min(A1),max(A1),length=50),yo=seq(min(A2),max(A2),length=50),duplicate='mean')
contour(INTERP2, levels=0.05, labels=' 95 % CI ', add = TRUE, col='#33333399', lwd=2, labcex=1.2)
abline(v=ENVonTCang[c(2,6)],col='#66ff00',lwd=4) # line for GMSST
axis(1,at=c(TCloc,270),labels=c('-E','F','A','I','E','-F','-A','-I','-E'),cex.axis=1.2)
axis(2,at=Lev,labels=p,cex.axis=1.2,las=1)
#axis(2,at=Lev,labels=c(1000,'',p[-c(1,2)]),cex.axis=1,las=1) # when using log scale p
axis(3,at =ENVonTCang,labels=c('-C','O','W','N','C','-O','-W','-N'),cex.axis=1.2)
axis(4, at=Lev,labels= rep('',length(Lev)),las=1)
mtext('TC Climate',side=1,line=3, cex=1.7)
mtext('Best Explanatory Environment',side=3,line=2.8, cex=1.7)
mtext('Pressure (hPa)',side=2,line=3.8, cex=1.6)
if (varnum==1) mtext( 'a',font=2,side=3,line=5.5, cex=3,adj=-0.2)
if (varnum==1) mtext( 'Moist Static Energy',font=1,side=3,line=5.5, cex=2,adj=-0.18)
if (varnum==2) mtext( 'b',font=2,side=3,line=5.5, cex=3,adj=-0.2)
if (varnum==2) mtext( 'Geopotential Height',font=1,side=3,line=5.5, cex=2,adj=-0.18)
box()
}
Fig. 3 Correlation maps of (a) moist static energy, and (b) geopotential height at each environmental variability direction (upper abscissa) in the tropics (30S-30N). TC climate variability which can be best explained by each environmental variability is noted in the lower abscissa. Green vertical lines indicate the variability direction of O and -O.
Observations indicate that global mean SST has increased by .30C over the past 30 years on average (Supplementary Fig. 9). Figure 4 shows how the increased global ocean warmth has influenced the traded-off between the intensity and the frequency of global TCs. Here, \( \sigma \) is the standard deviation of global intensity and frequency simultaneously, which suggests a trade-off rate as .21 ms\( ^{-1} \) per TC (Supplementary Fig. 12). Inversely, the rate results in 4.7 fewer TCs per meter per second increase in average intensity. The efficiency of intensity can be estimated on this \( \sigma \) level to estimate the amount of trade-off. The efficiency of intensity has increased by .74\( \sigma \) due to the global mean SST increase. In physical units the trade-off amounts to 1.3 ms\( ^{-1} \) at the expense of 6.1 TCs over the period.
In spite of the regionally strong response of the efficiency of intensity in the western North Pacific (see Fig. 2) to warming seas, the contribution of this region to the global result is small. This is because the intensity level in this basin is higher so a reduction in frequency implies a smaller contribution to the global intensity-frequency trade-off (Supplementary Fig. 5). Interestingly, the contribution in the North Atlantic is impressive. Despite the fact that the regional change of the efficiency of intensity is negative with global ocean warming (see Fig. 2), the contribution of this region to the global is positive and dominant. From the framework of statistical physics, the regional changes associated with global ocean warmth fix the global trade-off between intensity and frequency.
1.set option
Year = 1984:2012
# manually input trade_off_int_char & trade_off_frq_char in section # 5
2.read data
frqgl = read.table("Data/frqgl.txt", TRUE)[, 1]
intgl = read.table("Data/intgl.txt", TRUE)[, 1]
frqsp = read.table("Data/frqsp.txt", TRUE)[, 1]
intsp = read.table("Data/intsp.txt", TRUE)[, 1]
frqsi = read.table("Data/frqsi.txt", TRUE)[, 1]
intsi = read.table("Data/intsi.txt", TRUE)[, 1]
frqna = read.table("Data/frqna.txt", TRUE)[, 1]
intna = read.table("Data/intna.txt", TRUE)[, 1]
frqenp = read.table("Data/frqenp.txt", TRUE)[, 1]
intenp = read.table("Data/intenp.txt", TRUE)[, 1]
frqwnp = read.table("Data/frqwnp.txt", TRUE)[, 1]
intwnp = read.table("Data/intwnp.txt", TRUE)[, 1]
frqni = read.table("Data/frqni.txt", TRUE)[, 1]
intni = read.table("Data/intni.txt", TRUE)[, 1]
gmsst = read.table("Data/gmsst_12m_1854to2013.dat", TRUE)[(Year[1] - 1853):(rev(Year)[1] -
1853), 1]
3.GMSST change over 30 years (from linear trend)
GMSST_change = 30 * lm(gmsst ~ Year)$coeff[2]
4.Global E change by GMSST change
eintgl = (scale(intgl) - scale(frqgl))/sqrt(2)
E = eintgl
model = lm(E ~ gmsst)
d.f = length(Year) - 2
slp = summary(model)$coefficients[2, 1] #slope
sde = summary(model)$coefficients[2, 2] #standard error of the slope
confint95 = sde * qt(0.975, df = d.f) #confidence interval=sderror* tvalue(at 95% C.I with d.f)
eint_change = GMSST_change * slp
eint_changeU = GMSST_change * (slp + confint95) # 95% CI (upper limit)
eint_changeL = GMSST_change * (slp - confint95) # 95% CI (lower limit)
5.Regional contributions to global E change
# separation of regional portions
# frqgl=frqni + frqwnp + frqenp + frqna + frqsi + frqsp # global FRQ
# intgl=intni2 + intwnp2 + intenp2 + intna2 + intsi2 +intsp2 # global INT
intsp2=intsp*frqsp/frqgl # SP portion in global INT
intsi2= intsi*frqsi/frqgl # SI portion in global INT
intna2=intna*frqna/frqgl # NAportion in global INT
intenp2= intenp*frqenp/frqgl # ENP portion in global INT
intwnp2= intwnp*frqwnp/frqgl # WNP portion in global INT
intni2 = intni*frqni/frqgl # NI portion in global INT
#global EINT constitution
#eintgl=(scale(intgl) - scale(frqgl))/sqrt(2) #global EINT
# =intgl/sd(intgl)/sqrt(2)-mean(intgl)/sd(intgl)/sqrt(2)-frqgl/sd(frqgl)/sqrt(2)+mean(frqgl)/sd(frqgl)/sqrt(2)
# = (intni2 + intwnp2+ intenp2+ intna2 + intsi2 + intsp2)/sd(intgl)/sqrt(2)
# -(frqni + frqwnp + frqenp + frqna + frqsi + frqsp)/sd(frqgl)/sqrt(2)
# +mean(frqgl)/sd(frqgl)/sqrt(2) -mean(intgl)/sd(intgl)/sqrt(2)
# = P6 + P5 + P4 + P3 + P2 + P1 + Const
P1 = intsp2/sd(intgl)/sqrt(2) - frqsp/sd(frqgl)/sqrt(2) # SP variant part
P2 = intsi2/sd(intgl)/sqrt(2) - frqsi/sd(frqgl)/sqrt(2) # SI variant part
P3 = intna2/sd(intgl)/sqrt(2) - frqna/sd(frqgl)/sqrt(2) # NA variant part
P4 = intenp2/sd(intgl)/sqrt(2) - frqenp/sd(frqgl)/sqrt(2) # ENP variant part
P5 = intwnp2/sd(intgl)/sqrt(2) - frqwnp/sd(frqgl)/sqrt(2) # WNP variant part
P6 = intni2/sd(intgl)/sqrt(2) - frqni/sd(frqgl)/sqrt(2) # NI variant Part
#reginoal portions of the global EINT change
P_change=NULL
P_changeL=NULL ; P_changeU=NULL # including Lower & Upper 95 conf intervals
for (i in 1:6)
{ if (i==1) P=P1; if (i==2) P=P2; if (i==3) P=P3
if (i==4) P=P4; if (i==5) P=P5; if (i==6) P=P6
model=lm(P~gmsst)
d.f=length(Year)-2
slp=summary(model)$coefficients[2,1] #slope
sde=summary(model)$coefficients[2,2] #standard error of the slope
confint95=sde*qt(0.975,df=d.f) #confidence interval=sderror* tvalue(at 95% C.I with d.f)
P_change = c(P_change, GMSST_change*slp)
P_changeU= c(P_changeU, GMSST_change*(slp+confint95))
P_changeL= c(P_changeL, GMSST_change*(slp-confint95))
}
# confirm the process (check if eint_change and the sum of regional changes are the same)
eint_change[1] ; sum(P_change)
## [1] 1.044
# convert all values into reference sigma (standard deviation) scale
# (reference sigma: sigma of frequency or intensity)
eint_change2 =eint_change/sqrt(2)
eint_changeL2 =eint_changeL/sqrt(2)
eint_changeU2 =eint_changeU/sqrt(2)
P_change2 =P_change/sqrt(2)
P_changeL2 =P_changeL/sqrt(2)
P_changeU2 =P_changeU/sqrt(2)
# trade-off amounts of frequency and intensity in EINT changes
trade_off_int = # calculate trade-off amount of int # kt*0.514444--> m/s
0.514444*eint_change2*c(sd(intgl),sd(intsp),sd(intsi),sd(intna),sd(intenp),sd(intwnp),sd(intni))
trade_off_frq = # calculate trade-off amount of frq
eint_change2*c(sd(frqgl),sd(frqsp),sd(frqsi),sd(frqna),sd(frqenp),sd(frqwnp),sd(frqni))
# write in character format for plotting
# (if values are used, for example, 3.0 can be plotted as just 3.)
round(trade_off_int,1)
## [1] 1.3 4.4 3.0 3.2 3.8 2.5 5.2
round(trade_off_frq,1)
## [1] 6.1 2.7 2.3 3.7 3.4 3.4 1.5
trade_off_sign = c('pos', 'pos', 'pos', 'neg', 'pos', 'pos', 'pos')
trade_off_int_char= c('1.3', '1.8', '1.4', '2.6', '2.0', '3.7', '1.6')
trade_off_frq_char=c( '6.1', '1.1', '1.1', '3.0', '1.8', '5.0', '0.4')
6.plot
par(mfrow=c(1,1), mar=c(5,1,1,1))
# global
vect = as.vector(c(eint_change2[1],P_change2))
vectL = as.vector(c(eint_changeL2[1],P_changeL2))
vectU = as.vector(c(eint_changeU2[1],P_changeU2))
collist=c('#333333','#cccc66','#cc9966','#cc6666','#cc9999','#cc99cc','#ccccff')
hh=barplot(vect,border=NA,col=collist,horiz=T,axes=FALSE
,xlim=c(-2,2),ylim=c(0,10))
## Warning: the condition has length > 1 and only the first element will be used
## Warning: the condition has length > 1 and only the first element will be used
## Warning: the condition has length > 1 and only the first element will be used
segments(vectL,hh,vectU,hh,col='#666666') # 95% confidence interval
points(vect,hh,pch=16,cex=1.5,col='#666666')
text(vect[1]+0.55 ,hh[1]+0.37, bquote(.(trade_off_int_char[1])*" "*m*" "*s^-1),cex=1.4,col='#993333')
text(vect[1]+0.53 ,hh[1]-0.29 , bquote(per*" "*.(trade_off_frq_char[1])*" "*TCs),cex=1.4,col='#333366')
for (i in 1:7) # add regional portions in sigma
{if (vect[i]>=0) # positive contribution
{text(vect[i]+0.15,hh[i]+0.32, bquote(.(round(vect[i],2))*sigma),cex=1,col='#993333') }
if (vect[i]<0) # negative contribution
{text(vect[i]-0.15,hh[i]+0.32, bquote(.(round(vect[i],2))*sigma),cex=1,col='#333366')}
}
nameset=c('South Pacific','South Indian','North Atlantic'
,'eastern North Pacific','western North Pacific','North Indian')
text( rep(-1.3,6), hh[2:7], nameset,col=1,adj=c(0,-1.0),cex=1.2)
text(-1.3, hh[1], 'global net',col=1,adj=c(0,-1.0),cex=1.4,font=2)
axis(1, at = seq(-1,1.5, 0.5),labels=expression(-1~sigma,-0.5~sigma,0,0.5~sigma,1.0~sigma,1.5~sigma),cex.axis=1.2)
mtext(expression(paste('Reference ', sigma)), side=1, line=3.3, cex=1.3)
Fig. 4 Influence of global ocean warming on the trade-off between intensity and frequency of global tropical cyclones. 95% confidence intervals are shown in gray lines. Reference \( \sigma \) is the standard deviation of global intensity and frequency at the same time. Regional contributions are shown with different colors. Red and blue characters represent positive and negative changes, respectively. Global trade-off amount is seen to have increased by .74\( \sigma \) over the past 30 years, which is comparable to 1.3 ms\( ^{-1} \) increase of intensity at the expense of 6.1 TCs.
This study clarifies the connection between global TC climate and global warming. Instead of trend analysis, annual global mean SST is used to directly examine annual variations. This approach avoids the data quality problem occurring in a trend analysis when the observation technique advances over time. A continuous variability space is applied that provides a broader and deeper understanding of the relationships not accessible from a restricted focus on frequency and intensity alone.
Results show that global ocean warmth best explains the variability direction associated with the efficiency of intensity. Subsequently, global ocean warmth is confirmed to bring about convectively more unstable environment with more moisture in the lower troposphere, but a stronger high pressure aloft. This inhibits TC occurrences but provokes bursting intensities. Through the process, global warming over the ocean is seen to contribute to the efficiency of intensity. Intensity increases at the expense of frequency, hence no change in overall activity. These results have important implications for the theory of climate. The methodology can be applied to TCs generated by dynamical climate models to better understand the meaning of TC climatology in future scenarios.
Environmental variables come from the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP) reanalysis (www.esrl.noaa.gov/psd/data/gridded). To avoid the wind-conversion problem, global TC data set is organized by 1-minute wind observations from Joint Typhoon Warning Center (JTWC, www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks) and NOAA National Hurricane Center (NHC, www.nhc.noaa.gov/data). For the study, 29-year observations (1984~2012) are used. 1984 is the year when enhanced infrared (IR) satellite image analysis was developed, and also the year when Japan Meteorological Agency (JMA) employed objective Dvorak's satellite image analysis technique to the operation procedure (see www.wmo.int/pages/prog/www/tcp/documents/JMAoperationalTCanalysis.pdf), which is thought to have improved the climate consensus between JTWC and JMA observations. Considering the large number of western North Pacific TCs, the period of this study is also expected to increase the reliability of the global result.