Understanding the sensitivity of Limiting Intensity to SST

Sarah Strazzo, James B. Elsner, et al.

Question: Why are most GCMs unable to capture the observed sensitivity of limiting intensity (LI) to SST?

Is there a strong relationship between model horizontal resolution ad sensitivity?
Are modeled tropopause temperatures too warm?
Are GCM-generated TCs not intensified via the same physical processes as observed TCs? (i.e. are GCM-generated TCs even “TCs”?)

The plan:

Calculate sensitivity for as many GCMs as possible. Plot relationship between sensitivity and resolution.
Examine the spatial and temporal structure of tropoause temperatures for reanalysis (CFSR, ERA, MERRA) and each of the models (what to do about storage? I imagine these files will take up a fair amount of space…).
Use a multiple regression of LI onto SST and tropopause temperature for obs/reanalysis and for each of the models to try to understand the third one.
Maybe do a few case studies (?)

Step 1: Make the required packages available and source any necessary R scripts.

setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")

require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
require(plyr)

source("ohcFunctions.R")
source("selectData.R")
source("ismevplus.R")

Step 2: Acquire the SST data (using the function “reformat.SST” from ohcFunctions.R).

begin = 1981
end = 2009
ll = "+proj=longlat +datum=WGS84"
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"

nc = open.ncdf("C:/Users/Sarah/Dropbox/SpatialExtremes/sst.mnmean.nc")
avgsst = reformat.SST(nc, begin, end, 8, 9, 10)
rm(nc)

coordinates(avgsst) = c("lon", "lat")
proj4string(avgsst) = CRS(ll)
SST.sdf = spTransform(avgsst, CRS(lcc))

Step 3: Load and format MERRA 100 hPa temperature data.

As in Seidel et al. (2001, JGR) and Emanuel et al. (2013, J. Climate), 100 hPa temperatures are used as an estimate of the tropical tropopause. The data come from NASA's modern era reanalysis and were downloaded on 17 January 2014 from http://disc.sci.gsfc.nasa.gov/daac-bin/FTPSubset.pl?LOOKUPID_List=MAIMCPASM (analyzed state).

load("C:/Users/Sarah/Dropbox/SpatialExtremes/MERRAbyYear.RData")

years = avgtt[, 3:31]
tt = rowMeans(years)
lonlat = avgtt[, 1:2]
ASOtt = cbind(lonlat, tt)
colnames(ASOtt) = c("lon", "lat", "tt")
avgtt = data.frame(ASOtt)

coordinates(avgtt) = c("lon", "lat")
proj4string(avgtt) = CRS(ll)
TT.sdf = spTransform(avgtt, CRS(lcc))

Step 3a: Make sure the 100 hPa temperatures look realistic

plot(avgtt$lat, avgtt$tt, xlab = "Latitude", ylab = "100 hPa Temp (K)")

plot of chunk Fig1

Figure 1: Plot of latitude vs. 100 hPa temperature. As expected, temperatures increase as latitude increases

Step 4: Load best-track data

load("C:/Users/Sarah/Dropbox/SpatialExtremes/best.use.RData")
Wind.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]
Wind.df = subset(Wind.df, Type == "*")
Wind.df$WmaxS = Wind.df$WmaxS * 0.5144
Wind.df$Wmax = Wind.df$Wmax * 0.5144
Wind.df$maguv = Wind.df$maguv * 0.5144
# Wind.df = subset(Wind.df, WmaxS >= 33)
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
length(unique(Wind.df$Sid))  # 334 TCs (1979-2005)

## [1] 334


Wind.sdf = Wind.df
coordinates(Wind.sdf) = c("lon", "lat")
proj4string(Wind.sdf) = CRS(ll)
Wind.sdf = spTransform(Wind.sdf, CRS(lcc))

Step 5: Tessellate the basin with hexagons

set.seed(3042)
bb = bbox(Wind.sdf) * 1.2
bb[2, 1] = bbox(Wind.sdf)[2, 1] * 0.7
hpt = spsample(Wind.sdf, type = "hexagonal", n = 120, bb = bb, offset = c(0.5, 
    0.5))
hpg = HexPoints2SpatialPolygons(hpt)
np = length(hpg@polygons)
area = hpg@polygons[[1]]@area * 10^-6  # square km
area  # 837382.8

## [1] 837383

Step 6: Overlay the best-track, sst, and tt data onto the hexagons and calculate per hexagon statistics

hexidOBS = over(x = Wind.sdf, y = hpg)
Wind.sdf$hexid = hexidOBS

OBSstats = ddply(Wind.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid), 
    WmaxS = max(WmaxS))
OBSstats2 = ddply(OBSstats, .(hexid), summarize, countOBS = length(Sid), WmaxOBS = max(WmaxS))

hexData = subset(OBSstats2, countOBS > 20)
hexData = hexData[hexData$hexid != 49, ]  # (over land)
hexData = hexData[hexData$hexid != 50, ]  # (over land)
hexData = hexData[hexData$hexid != 25, ]  # (over land)
row.names(hexData) = paste("ID", hexData$hexid, sep = "")
hspdf = SpatialPolygonsDataFrame(hpg[row.names(hexData)], hexData)

avgsst = over(x = hpg[row.names(hexData)], y = SST.sdf, fn = mean)
hspdf$avgsst = avgsst$avgsst

avgtt = over(x = hpg[row.names(hexData)], y = TT.sdf, fn = mean)
hspdf$tt = avgtt$tt

Step 7: Set up map borders

cl = map("worldHires", xlim = c(-120, 20), ylim = c(-10, 70), plot = TRUE, resolution = 0.7)
clp = map2SpatialLines(cl, proj4string = CRS(ll))
clp = spTransform(clp, CRS(lcc))
l1 = list("sp.lines", clp, col = "gray")

cr1 = brewer.pal(6, "Greens")
rng1 = seq(15, 75, 10)
p0 = spplot(hspdf, "countOBS", col = "white", col.regions = cr1, sp.layout = list(l1), 
    at = rng1, colorkey = list(space = "bottom", labels = paste(rng1)), sub = expression("Observed TC Frequency (1981-2009)"))

cr2 = brewer.pal(6, "Oranges")
rng2 = seq(25, 85, 10)
p1 = spplot(hspdf, "WmaxOBS", col = "white", col.regions = cr2, sp.layout = list(l1), 
    at = rng2, colorkey = list(space = "bottom", labels = paste(rng2)), sub = expression("Per Hexagon Max Wind Speed (m/s)"))

cr3 = brewer.pal(7, "Reds")
rng3 = seq(9, 30, 3)
p2 = spplot(hspdf, "avgsst", col = "white", col.regions = cr3, sp.layout = list(l1), 
    at = rng3, colorkey = list(space = "bottom", labels = paste(rng3)), sub = expression("Mean ASO SST"))

cr4 = rev(brewer.pal(6, "Blues"))
rng4 = seq(194, 218, 4)
p3 = spplot(hspdf, "tt", col = "white", col.regions = cr4, sp.layout = list(l1), 
    at = rng4, colorkey = list(space = "bottom", labels = paste(rng4)), sub = expression("Mean 100 hPa Temperature"))

p0 = update(p0, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p1 = update(p1, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p2 = update(p2, main = textGrob("c", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p3 = update(p3, main = textGrob("d", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))

plot(p0, split = c(1, 1, 2, 2), more = TRUE)
plot(p1, split = c(2, 1, 2, 2), more = TRUE)
plot(p2, split = c(1, 2, 2, 2), more = TRUE)
plot(p3, split = c(2, 2, 2, 2), more = FALSE)

plot of chunk Fig2

Figure 2: TC Frequency, per hexagon maximum TC intensity, SST, and 100 hPa temperature maps

Step 8: Explore the data

data = hspdf@data
cor.test(data$avgsst, data$tt)

## 
##  Pearson's product-moment correlation
## 
## data:  data$avgsst and data$tt
## t = -12.36, df = 28, p-value = 7.361e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9612 -0.8359
## sample estimates:
##     cor 
## -0.9193


cor.test(data$avgsst, data$countOBS)

## 
##  Pearson's product-moment correlation
## 
## data:  data$avgsst and data$countOBS
## t = 1.781, df = 28, p-value = 0.0858
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04668  0.60922
## sample estimates:
##   cor 
## 0.319

cor.test(data$avgsst, data$WmaxOBS)

## 
##  Pearson's product-moment correlation
## 
## data:  data$avgsst and data$WmaxOBS
## t = 4.444, df = 28, p-value = 0.0001265
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3682 0.8147
## sample estimates:
##    cor 
## 0.6431


cor.test(data$tt, data$countOBS)

## 
##  Pearson's product-moment correlation
## 
## data:  data$tt and data$countOBS
## t = -0.6813, df = 28, p-value = 0.5013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4665  0.2438
## sample estimates:
##     cor 
## -0.1277

cor.test(data$tt, data$WmaxOBS)

## 
##  Pearson's product-moment correlation
## 
## data:  data$tt and data$WmaxOBS
## t = -3.669, df = 28, p-value = 0.001014
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7716 -0.2636
## sample estimates:
##     cor 
## -0.5698

A few notes:

ASO SST and ASO 100 hPa temperatures are strongly negatively correlated…
As with ASO SST, ASO 100 hPa temperatures are not very highly correlated to TC frequency
ASO 100 hPa temperatures are negatively correlated

xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("100 hPa temperature [", degree, "C]"))
ggplot(data[data$avgsst > 25, ], aes(x = avgsst, y = tt)) + geom_point() + geom_smooth(method = "lm") + 
    xlab(xlabel) + ylab(ylabel) + theme_bw()

plot of chunk Fig3

Figure 3: Per hexagon ASO SST vs. ASO 100 hPa temperatures. Only ASO SST > 25 C are shown

Step 9: Obtain estimates for limiting intensity using GPD/Poisson model.

oDat = OBSstats[OBSstats$hexid != 49, ]
oDat = oDat[oDat$hexid != 50, ]
oDat = oDat[oDat$hexid != 25, ]
oDat = split(oDat, oDat$hexid)
less = which(sapply(oDat, function(x) dim(x)[1] < 21))
oDat = oDat[-less]

m = length(begin:end)
lambda = numeric()
sigma = numeric()
xi = numeric()
thresh = 30
rate = numeric()
theo = numeric()

for (i in 1:length(oDat)) {
    W = oDat[[i]]$WmaxS
    rate[i] = length(W)/m
    model = exceed.pp.plot(W, npy = rate[i], threshold = thresh, yunits = "m/s", 
        addtext = "", plotrl = FALSE)
    lambda[i] = model$gpd.coef[1]
    sigma[i] = model$gpd.coef[2]
    xi[i] = model$gpd.coef[3]
    theo[i] = thresh - sigma[i]/xi[i]
}

Step 9a: Examine model parameters. Make sure everything looks reasonable.

hspdf$lambda = lambda
hspdf$sigma = sigma
hspdf$xi = xi
hspdf$theo = theo

xpo = sum(xi > 0)  # All models are asymptotic
xl1 = sum(xi < -1)  # 16/30 LIs = per hexagon WmaxS.  Not stellar...

cr1 = brewer.pal(7, "Greens")
rng1 = seq(0, 1.4, 0.2)
p0 = spplot(hspdf, "lambda", col = "white", col.regions = cr1, sp.layout = list(l1), 
    at = rng1, colorkey = list(space = "bottom", labels = paste(rng1)), sub = expression(paste("Parameter: ", 
        lambda)))

cr2 = brewer.pal(7, "Oranges")
rng2 = seq(0, 70, 10)
p1 = spplot(hspdf, "sigma", col = "white", col.regions = cr2, sp.layout = list(l1), 
    at = rng2, colorkey = list(space = "bottom", labels = paste(rng2)), sub = expression(paste("Paramter: ", 
        sigma)))

cr3 = brewer.pal(8, "Reds")
rng3 = seq(-4, 0, 0.5)
p2 = spplot(hspdf, "xi", col = "white", col.regions = cr3, sp.layout = list(l1), 
    at = rng3, colorkey = list(space = "bottom", labels = paste(rng3)), sub = expression(paste("Paramter: ", 
        xi)))

cr4 = brewer.pal(5, "Blues")
rng4 = seq(30, 80, 10)
p3 = spplot(hspdf, "theo", col = "white", col.regions = cr4, sp.layout = list(l1), 
    at = rng4, colorkey = list(space = "bottom", labels = paste(rng4)), sub = expression("LI"))

p0 = update(p0, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p1 = update(p1, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p2 = update(p2, main = textGrob("c", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p3 = update(p3, main = textGrob("d", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))

plot(p0, split = c(1, 1, 2, 2), more = TRUE)
plot(p1, split = c(2, 1, 2, 2), more = TRUE)
plot(p2, split = c(1, 2, 2, 2), more = TRUE)
plot(p3, split = c(2, 2, 2, 2), more = FALSE)

plot of chunk Fig4

Figure 4: Model parameters

Step 10: Try to calculate sensitiviy….

dat = as.data.frame(hspdf)
dat$northing = coordinates(hspdf)[, 2]
dat = dat[avgsst >= 25, ]


summary(lm(theo ~ avgsst, weights = countOBS, data = dat))

## 
## Call:
## lm(formula = theo ~ avgsst, data = dat, weights = countOBS)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -175.60  -25.24    1.31   30.03  100.87 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -155.60      53.21   -2.92  0.00786 ** 
## avgsst          7.83       1.90    4.12  0.00045 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57.9 on 22 degrees of freedom
## Multiple R-squared:  0.435,  Adjusted R-squared:  0.409 
## F-statistic: 16.9 on 1 and 22 DF,  p-value: 0.000454

summary(lm(theo ~ tt, weights = countOBS, data = dat))

## 
## Call:
## lm(formula = theo ~ tt, data = dat, weights = countOBS)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -190.88  -54.60    3.73   52.35  101.90 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  177.996    158.124    1.13     0.27
## tt            -0.576      0.794   -0.73     0.48
## 
## Residual standard error: 76.2 on 22 degrees of freedom
## Multiple R-squared:  0.0234, Adjusted R-squared:  -0.021 
## F-statistic: 0.526 on 1 and 22 DF,  p-value: 0.476

summary(lm(theo ~ avgsst + tt, weights = countOBS, data = dat))

## 
## Call:
## lm(formula = theo ~ avgsst + tt, data = dat, weights = countOBS)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -131.19  -28.86    1.57   19.02   72.67 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -747.936    189.913   -3.94  0.00075 ***
## avgsst        12.831      2.230    5.75    1e-05 ***
## tt             2.272      0.708    3.21  0.00422 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 48.6 on 21 degrees of freedom
## Multiple R-squared:  0.621,  Adjusted R-squared:  0.585 
## F-statistic: 17.2 on 2 and 21 DF,  p-value: 3.77e-05


# lmttSST = summary(lm(tt ~ avgsst, weights=countOBS, data=dat)) resTT =
# residuals(lmttSST) lmLIttRes = summary(lm(dat$theo ~ resTT,
# weights=dat$countOBS))

Per grid ASO 100 hPa temperatures are highly correlated to SST. Adding tropopause temperature as a covariate in the sensitivity model does not seem wise given the high correlation between the two variables.

There is no significant relationship between 100 hPa temperatures and limiting intensity. Perhaps this is because atmospheric temperatures are less homogeneous compared to SST? They change on shorter timescales and, overall, are more spatially local than SST. Therefore, averaging over three months and over a large area (size of a hexagon) might not be useful here. Is there a better way to examine their influence? What does this mean about the effect of ENSO on LI?

Perhaps I'll try using MDR tropospheric temperatures and quantile regression. How would that work if you accounted for ENSO, though? How much of the year-to-year variability in tropopause tempeartures is caused by ENSO?

But because trends in tropopause temps have come up on several occasions, I do want to quickly examine spatial patterns in these trends over the basin.

How have tropopause temperatures changed in time?

Step 1: Re-open MERRAbyYear data, assign a CRS, and overlay on the grid.

load("C:/Users/Sarah/Dropbox/SpatialExtremes/MERRAbyYear.RData")
byYear = avgtt
colnames(byYear) = c("lon", "lat", "t1981", "t1982", "t1983", "t1984", "t1985", 
    "t1986", "t1987", "t1988", "t1989", "t1990", "t1991", "t1992", "t1993", 
    "t1994", "t1995", "t1996", "t1997", "t1998", "t1999", "t2000", "t2001", 
    "t2002", "t2003", "t2004", "t2005", "t2006", "t2007", "t2008", "t2009")
coordinates(byYear) = c("lon", "lat")
proj4string(byYear) = CRS(ll)
TTy.sdf = spTransform(byYear, CRS(lcc))

tty = over(x = hpg[row.names(hexData)], y = TTy.sdf, fn = mean)

Step 2: Calculate the per grid trend using a linear regression of 100 hPa temperature onto year (Perhaps I should account for ENSO?).

years = 1981:2009
trends = matrix(0, ncol = 2, nrow = dim(tty)[1])
trends = data.frame(trends)

for (i in 1:dim(tty)[1]) {
    temps = tty[i, ]
    temps = t(temps)
    rownames(temps) = 1:29
    mdat = cbind(years, temps)
    mdat = data.frame(mdat)
    colnames(mdat) = c("years", "temps")
    mod = lm(temps ~ years, data = mdat)
    slope = coef(mod)[2]
    se = sqrt(diag(vcov(mod)))[2]
    trends[i, 1] = slope
    trends[i, 2] = se
}

colnames(trends) = c("slope", "se")
rownames(trends) = rownames(tty)
hspdf$trend = trends$slope
hspdf$se = trends$se

centers = coordinates(hspdf)
l2 = list("sp.text", centers, round(hspdf$se, 3), col = "red", cex = 1.1)

cr = rev(brewer.pal(8, "Blues"))
rng = seq(-0.1, -0.02, 0.01)
spplot(hspdf, "trend", col = "white", col.regions = cr, sp.layout = list(l1, 
    l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Change in 100 hPa Temp per Year"))

plot of chunk Fig5

Figure 5: Change in 100 hPa temperature per year (1981–2009) Red text indicates standard error

How does ENSO affect this relationship? (coming soon)

con = "http://hurricaneclimatology.squarespace.com/storage/chapter-5/SOI.txt"
SOI = read.table(con, header = TRUE)
SOI = subset(SOI, Year > 1980 & Year < 2010)
aso = SOI[, 8:10]
aso = rowMeans(aso)
soi = cbind(SOI$Year, aso)
colnames(soi) = c("Year", "ASO")