Using the current methodology, we average SSTs over both space and time. We take the area-averaged per hexagon ASO SST that is also temporally averaged over the 1979–2010 time period. Would it be better to use a more temporally local SST? Here's what I propose:
Justification: Why am I bothering with temporally local SSTs rather than ASO averaged SSTs over the entire 1980–2010 time period?
This allows me to calculate “relative SST,” which I think some folks might be interested in seeing. (Spoiler alert: Interestingly, the sensitivity is approximately the same for relative vs. “lcoal” SST)
This (potentially) allows me to examine not simply the “stationary” spatial relationship between SST and LI, but also incorporate some temporal aspect, although I admittedly haven't decided how to do this yet.
Note: It seems possible that in some cases, adjacent hexagons will have maximum wind speeds from the same TC (i.e. We can't assume independence). This seems like a problem….
Note for Dr. Elsner: Be forewarned that there are a lot of “for” loops in the code that follows! I'll remedy this later if any of it ends up being a worthwhile pursuit (or at the very least, turn the loop chunks into separate functions to clean things up a bit)!
setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")
require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
require(plyr)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
source("C:/Users/Sarah/Dropbox/SpatialExtremes/selectData.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/ismevplus.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Kerrycanes/sensPoly.R")
source("ohcFunctions.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensitivityUncertainty.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensTrends.R")
Load best track data, subset on tropical cycloens, convert wind speeds to meters per second, and estimate the rotational wind speed by subtracting 60% of the forward motion (maguv)
load("C:/Users/Sarah/Dropbox/Research/Data/best.use2.RData")
### Ruh roh...what's up with the dates not matching the TCs?
begin = 1979
end = 2010
Wind.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]
Wind.df = subset(Wind.df, Type == "TS" | Type == "HU")
Wind.df$WmaxS = Wind.df$WmaxS * 0.5144
quantile(Wind.df$WmaxS, prob = seq(0.6, 0.8, 0.01))
## 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71%
## 33.31 33.44 33.57 33.81 34.21 34.67 35.10 35.60 35.98 36.29 36.80 37.33
## 72% 73% 74% 75% 76% 77% 78% 79% 80%
## 37.87 38.36 38.67 39.11 39.80 40.52 41.09 41.65 42.49
Wind.df$maguv = Wind.df$maguv * 0.5144
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
length(unique(Wind.df$Sid)) # 372
## [1] 372
Load SST data and reformat it (using reformat.SST.temp) to generate a data frame of monthly SSTs (JJASO) for the years 1979–2010. Use relative.SST to generate a data frame of monthly tropical average (30N to 30S) SSTs for the same period. Assign a lambert conformal conic projection to each of the TC and SST data frames
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")
mSST = reformat.SST.temp(nc, 1979, 2010) ## (monthly SST)
tSST = relative.SST(nc, 1979, 2010) ## (tropical SST)
tSST = data.frame(tSST)
rm(nc)
coordinates(mSST) = c("lon", "lat")
proj4string(mSST) = CRS(ll)
mSST.sdf = spTransform(mSST, CRS(lcc))
Wind.sdf = Wind.df
coordinates(Wind.sdf) = c("lon", "lat")
proj4string(Wind.sdf) = CRS(ll)
Wind.sdf = spTransform(Wind.sdf, CRS(lcc))
Define a set of equal area (~543 515 km2) hexagons that completely cover the set of TC tracks
set.seed(1025)
bb = matrix(0, nrow = 2, ncol = 2)
bb[2, 1] = bbox(Wind.sdf)[2, 1] * 0.8
bb[1, 1] = bbox(Wind.sdf)[1, 1] * 1.35
bb[1, 2] = bbox(Wind.sdf)[1, 2] * 1.08
bb[2, 2] = bbox(Wind.sdf)[2, 2] * 1.08
hpt = spsample(Wind.sdf, type = "hexagonal", n = 150, bb = bb)
hpg = HexPoints2SpatialPolygons(hpt)
# plot(hpg) plot(Wind.sdf, add=TRUE, pch='.', col='red')
np = length(hpg@polygons) # 122
area = hpg@polygons[[1]]@area * 10^-6
# 543514.6 km^2
Overlay the track data onto the hexagon data and exract the entire record associated with the per hexagon maximum wind speed
hexid = over(x = Wind.sdf, y = hpg)
Wind.sdf$hexid = hexid
windSplit = split(Wind.sdf@data, Wind.sdf$hexid)
maxWinds = matrix(0, nrow = length(windSplit), ncol = dim(Wind.sdf@data)[2])
maxWinds = data.frame(maxWinds)
colnames(maxWinds) = colnames(Wind.sdf@data)
r = 1
for (i in 1:length(windSplit)) {
dat = windSplit[[i]]
m = max(dat$WmaxS)
m2 = subset(dat, WmaxS == m)
maxWinds[r, ] = m2
r = r + 1
}
length(table(maxWinds$name)) # 41 storms for 82 grids
## [1] 41
plot(table(maxWinds$name))
Calculate the per hexagon TC count and subset on only those hexagons with at least 15 TCs
stats = ddply(Wind.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid), WmaxS = max(WmaxS))
stats2 = ddply(stats, .(hexid), summarize, countOBS = length(Sid))
maxWinds$count = stats2$countOBS
hexData = subset(maxWinds, count >= 15)
hexData = subset(hexData, hexid != 72) # Why did I get rid of this one??
row.names(hexData) = paste("ID", hexData$hexid, sep = "")
mSST = over(x = hpg[row.names(hexData)], y = mSST.sdf, fn = mean)
Find the per hexagon monthly average SST corresponding to the month/year of each per hexagon maximum wind speed
localsst = numeric()
for (i in 1:dim(hexData)[1]) {
month = hexData$Mo[i]
if (month < 10) {
date = paste("Y", hexData$Yr[i], "M", formatC(month, 1, flag = "0"),
sep = "")
} else {
date = paste("Y", hexData$Yr[i], "M", formatC(month), sep = "")
}
s = mSST[i, colnames(mSST) == date]
localsst[i] = s
}
localsst = data.frame(localsst)
rownames(localsst) = rownames(hexData)
hexData$localsst = localsst
Find the per hexagon monthly average SST corresponding to the month prior to each recorded per hexagon maximum wind speed (same year)
localsst2 = numeric()
for (i in 1:dim(hexData)[1]) {
month = hexData$Mo[i] - 1
if (month < 10) {
date = paste("Y", hexData$Yr[i], "M", formatC(month, 1, flag = "0"),
sep = "")
} else {
date = paste("Y", hexData$Yr[i], "M", formatC(month), sep = "")
}
s = mSST[i, colnames(mSST) == date]
localsst2[i] = s
}
localsst2 = data.frame(localsst2)
rownames(localsst2) = rownames(hexData)
hexData$localsst2 = localsst2
Find the tropical average SST (30N to 30S) for the month of each recorded per hexagon maximum intensity
tropSST = numeric()
for (i in 1:dim(hexData)[1]) {
month = hexData$Mo[i]
if (month < 10) {
date = paste("Y", hexData$Yr[i], "M", formatC(month, 1, flag = "0"),
sep = "")
} else {
date = paste("Y", hexData$Yr[i], "M", formatC(month), sep = "")
}
s2 = tSST[rownames(tSST) == date, 1]
tropSST[i] = s2
}
tropSST = data.frame(tropSST)
rownames(tropSST) = rownames(hexData)
hexData$tropSST = tropSST
Find the tropical average SST for the month prior to each per hexagon maximum intensity (same year)
tropSST2 = numeric()
for (i in 1:dim(hexData)[1]) {
month = hexData$Mo[i] - 1
if (month < 10) {
date = paste("Y", hexData$Yr[i], "M", formatC(month, 1, flag = "0"),
sep = "")
} else {
date = paste("Y", hexData$Yr[i], "M", formatC(month), sep = "")
}
s2 = tSST[rownames(tSST) == date, 1]
tropSST2[i] = s2
}
tropSST2 = data.frame(tropSST2)
rownames(tropSST2) = rownames(hexData)
hexData$tropSST2 = tropSST2
Calcuate the average ASO tropical SST (30N to 30S) corresponding to the year of each recorded per hexagon maximum intensity
ASOtropSST = numeric()
aug = 8
sep = 9
oct = 10
for (i in 1:dim(hexData)[1]) {
m1 = paste("Y", hexData$Yr[i], "M", formatC(aug, 1, flag = "0"), sep = "")
m2 = paste("Y", hexData$Yr[i], "M", formatC(sep, 1, flag = "0"), sep = "")
m3 = paste("Y", hexData$Yr[i], "M", formatC(oct), sep = "")
A = tSST[rownames(tSST) == m1, 1]
S = tSST[rownames(tSST) == m2, 1]
O = tSST[rownames(tSST) == m3, 1]
ASOtropSST[i] = (A + S + O)/3
}
ASOtropSST = data.frame(ASOtropSST)
rownames(ASOtropSST) = rownames(hexData)
hexData$ASOtropSST = ASOtropSST
Find the per hexagon average ASO SST corresponding to the year of each per hexagon maximum intensity
ASOsst = numeric()
aug = 8
sep = 9
oct = 10
for (i in 1:dim(hexData)[1]) {
m1 = paste("Y", hexData$Yr[i], "M", formatC(aug, 1, flag = "0"), sep = "")
m2 = paste("Y", hexData$Yr[i], "M", formatC(sep, 1, flag = "0"), sep = "")
m3 = paste("Y", hexData$Yr[i], "M", formatC(oct), sep = "")
A = mSST[i, colnames(mSST) == m1]
S = mSST[i, colnames(mSST) == m2]
O = mSST[i, colnames(mSST) == m3]
ASOsst[i] = (A + S + O)/3
}
ASOsst = data.frame(ASOsst)
rownames(ASOsst) = rownames(hexData)
hexData$ASOsst = ASOsst
Same as above, but for average JAS SST
JASsst = numeric()
jul = 7
aug = 8
sep = 9
for (i in 1:dim(hexData)[1]) {
m1 = paste("Y", hexData$Yr[i], "M", formatC(jul, 1, flag = "0"), sep = "")
m2 = paste("Y", hexData$Yr[i], "M", formatC(aug, 1, flag = "0"), sep = "")
m3 = paste("Y", hexData$Yr[i], "M", formatC(sep, 1, flag = "0"), sep = "")
J = mSST[i, colnames(mSST) == m1]
A = mSST[i, colnames(mSST) == m2]
S = mSST[i, colnames(mSST) == m3]
JASsst[i] = (J + A + S)/3
}
JASsst = data.frame(JASsst)
rownames(JASsst) = rownames(hexData)
hexData$JASsst = JASsst
Same as above, but for average JJA SST
JJAsst = numeric()
jun = 6
jul = 7
aug = 8
for (i in 1:dim(hexData)[1]) {
m1 = paste("Y", hexData$Yr[i], "M", formatC(jun, 1, flag = "0"), sep = "")
m2 = paste("Y", hexData$Yr[i], "M", formatC(jul, 1, flag = "0"), sep = "")
m3 = paste("Y", hexData$Yr[i], "M", formatC(aug, 1, flag = "0"), sep = "")
J1 = mSST[i, colnames(mSST) == m1]
J2 = mSST[i, colnames(mSST) == m2]
A = mSST[i, colnames(mSST) == m3]
JJAsst[i] = (J1 + J2 + A)/3
}
JJAsst = data.frame(JJAsst)
rownames(JJAsst) = rownames(hexData)
hexData$JJAsst = JJAsst
Calculate the per hexagon average ASO SST over the entire 1979–2010 period
allSST = mSST[, 65:160]
avgsst = rowMeans(allSST)
avgsst = data.frame(avgsst)
hexData$avgsst = avgsst
Calculate relative SST (local-tropical mean) for the month of and month prior to the recorded per hexagon maximum intensity
relsst = hexData$localsst - hexData$tropSST
rownames(relsst) = rownames(hexData)
colnames(relsst) = c("relsst")
hexData$relsst = relsst
relsst2 = hexData$localsst2 - hexData$tropSST2
rownames(relsst2) = rownames(hexData)
colnames(relsst2) = c("relsst2")
hexData$relsst2 = relsst2
ASOrelSST = hexData$ASOsst - hexData$ASOtropSST
rownames(ASOrelSST) = rownames(hexData)
colnames(ASOrelSST) = c("ASOrelSST")
hexData$ASOrelSST = ASOrelSST
Create a spatial polygons data frame
hspdf = SpatialPolygonsDataFrame(hpg[row.names(hexData)], hexData)
Combine the intensity, count and assorted SST data
hd = cbind(hspdf$WmaxS, hspdf$localsst, hspdf$relsst, hspdf$localsst2, hspdf$relsst2,
hspdf$avgsst, hspdf$count, hspdf$ASOsst, hspdf$JASsst, hspdf$JJAsst, hspdf$ASOrelSST,
hspdf$Yr)
colnames(hd) = c("WmaxS", "localsst", "relsst", "localPre", "relPre", "avgsst",
"count", "ASO", "JAS", "JJA", "ASOrel", "Yr")
Calculate correlation between per hexagon maximum intensity and assorted SSTs
cor.test(hd$WmaxS, hd$localsst)
# 0.674 (0.482, 0.804)
cor.test(hd$WmaxS, hd$relsst)
# 0.675 (0.484, 0.805)
cor.test(hd$WmaxS, hd$localPre)
# 0.698 (0.516, 0.819)
cor.test(hd$WmaxS, hd$relPre)
# 0.698 (0.513, 0.818)
cor.test(hd$WmaxS, hd$avgsst)
# 0.686 (0.499, 0.812)
cor.test(hd$WmaxS, hd$ASO)
# 0.681 (0.492, 0.809)
cor.test(hd$WmaxS, hd$ASOrel)
# 0.682 (0.493, 0.809)
cor.test(hd$WmaxS, hd$JAS)
# 0.694 (0.512, 0.817)
cor.test(hd$WmaxS, hd$JJA)
# 0.698 (0.516, 0.820)
Regress per hexagon maximum intensity onto per hexagon average ASO SST for the 1979–2010 period
lmAS = lm(WmaxS ~ avgsst, weights = count, data = hd[hd$avgsst >= 25, ])
summary(lmAS)
# 8.88 +/- 1.42 R2 = 0.529 AIC = 258.79
Note: This number stays fairly stable for SST thresholds ranging from 24-27C
Regress per hexagon maximum intensity onto per hexagon temporally local SST and relative SST
lmLS = lm(WmaxS ~ localsst, weights = count, data = hd[hd$localsst >= 25, ])
summary(lmLS)
# 6.49 +/- 1.16 R2 = .4577 AIC = 280.54
lmRS = lm(WmaxS ~ relsst, weights = count, data = hd[hd$localsst >= 25, ])
summary(lmRS)
# 6.23 +/- 1.13 R2 = .4494 AIC = 281.13
Why does average ASO SST (averaged over the entire 30-year period) fit the per hexagon max intensity better than local SST or relative SST? Could it be the mixing effect? Does the average SST for that hexagon during the month of the recorded max intensity decrease due to TC-induced mixing? I can test this by using local SSTs from the month prior to the date of the per hexagon maximum intensity.
Regress per hexagon maximum intensity onto per hexagon SST from the month prior to the per hexagon maximum intensity
lmLS2 = lm(WmaxS ~ localPre, weights = count, data = hd[hd$localPre >= 25, ])
summary(lmLS2)
# 6.74 +/- 1.32 R2 = .4217 AIC = 272.98
lmRS2 = lm(WmaxS ~ relPre, weights = count, data = hd[hd$localPre >= 25, ])
summary(lmRS2)
# 7.01 +/- 1.36 R2 = .4240 AIC = 272.83
Regress per hexagon maximum intensity onto per hexagon local ASO SST (from the year of the recorded per hexagon max intensity)
lmASO = lm(WmaxS ~ ASO, weights = count, data = hd[hd$ASO >= 25, ])
summary(lmASO)
# 10.29 +/- 1.46 R2 = .5863 AIC = 253.95
lmASOrel = lm(WmaxS ~ ASOrel, weights = count, data = hd[hd$ASO >= 25, ])
summary(lmASOrel)
# 10.04 +/- 1.44 R2 = .5803 AIC = 254.47
Plot the sensitivity of per hexagon maximum intensity to SST
C = hspdf@data
C = subset(C, avgsst > 25)
C1 = cbind(C$hexid, C$count, C$avgsst, C$WmaxS)
colnames(C1) = c("hexid", "count", "sst", "Vmax")
C1 = data.frame(C1)
C2 = cbind(C$hexid, C$count, C$localsst, C$WmaxS)
colnames(C2) = c("hexid", "count", "sst", "Vmax")
C2 = data.frame(C2)
C3 = cbind(C$hexid, C$count, C$localsst2, C$WmaxS)
colnames(C3) = c("hexid", "count", "sst", "Vmax")
C3 = data.frame(C3)
C4 = cbind(C$hexid, C$count, C$ASOsst, C$WmaxS)
colnames(C4) = c("hexid", "count", "sst", "Vmax")
C4 = data.frame(C4)
pdat = rbind(C1, C2, C3, C4)
pdat$type = factor(c(rep("ASO 1979-2010", dim(C)[1]), rep("Month Local", dim(C)[1]),
rep("Prior Month Local", dim(C)[1]), rep("ASO Local", dim(C)[1])))
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Per grid max intensity [m ", s^-1, "]"))
p25 = ggplot(pdat, aes(x = sst, y = Vmax)) + geom_point() + facet_grid(~type)
p25 = p25 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p25 = p25 + geom_smooth(method = "lm", aes(weight = count))
p25
There is very little difference between the four SST options; however, local ASO seems to be the closest fit (i.e. highest R2, lowest AIC, also the highest sensitivity). This also allows us to calculate the sensitivity of per hexagon maximum intensity to relative SST. Note that I also examined the sensitivity of per hexagon maximum intensity to local JAS and JJA SST but did not see as high a sensitivity or low an AIC.
To summarize:
ASO 1979–2010
Local (same month and year)
Relative (same month and year)
Local (prior month, same year)
Relative (prior month, same year)
Local ASO (same year)
Relative ASO (same year)
A few observations, thoughts, questins:
Why do temporally averaged SSTs explain more of the variance in max winds than do temporally local(ish) SSTs? (mixing? What happens if I look at average monthly SST as a function of the number of TCs that occured in that hexagon during that month? Perhaps that would shed some light on the lack of a relationship between maximum wind speed and temporally/spatially local SST.)
What if we examine temporal trends per grid (Note: Did this, saw nothing)? And sensitivity per grid? That might be interesting (although I doubt that the signal is strong…).
Why do we select 25C? Why not 26.5 or 26 or 24? (checked this: the sensitivity is fairly stable from about 24C to 27C. Above 27C there are too few data to estimate sensitivity. Below 24C, the relationship weakens.).
Top: Per hexagon maximum intensity over the 1979–2010 time period (red text indicates the name of the TC). Bottom: Average ASO SST for the month during which the maximum TC intensity was observed.
Of course, there are problems associated with this method:
Hexagons aren't statistically independent when you use the per hexagon maximum intensity (Although I suppose they weren't technically independent when we estimated LI using the GPD-Poisson distribution either). A per hexagon maximum wind speed may come from the same TC for multiple hexagons (see above figure). How does this affect the analysis? I can allow each TC to accont for the per hexagon max intensity only once, but is that an appropriate way to deal with this issue?? Besides TC intensity, per hexagon SST isn't truly independent either. Is there some way to account for this autocorrelation? Spatial regression perhaps? (Look at old code for help with this)
There's also the TC-induced mixing issue. As an intense TC moves over the area covered by a given hexagon, it will induce vertical mixing that effectively reduces the local SST for as long as a month following TC passage. I think this may be what we see when we use local SST in place of average ASO SST. The sensitivity is lower, the R2 is smaller, and the AIC is higher when we use local SST (“local” meaning average monthly SST from the month/year during which the per hexagon maximum intensity was recorded). If I use the month prior to the maximum intensity, the relationship improves but is still not as strong as average ASO SST. Note that average ASO SST from the year of the per hexagon maximum intensity explains more variance in intensity than does the 1979–2010 ASO average. This makes sense, thankfully. I assume ASO is better than local b/c we are averaging out (in a sense) the mixing effects.
Note: I started working on a “per hexagon sensitivity of per TC max intensity to SST,” but have not yet finished (plus, there are some independence issues that should probably be resolved first).
The maximum potential intensity of TCs is thought to be driven primarily by thermodynamic processes (e.g., Emanuel 1986; Holland 1997). Of particular importance is the enthalpy gradient between the sea surface and the atmosphere, which is partially dependent on SST. In light of this thermodynamically driven model for TC maximum potential intensity (MPI), most climate studies focus on the influence of local SST on TCs. Climate modeling groups adjust SSTs and examine the effects of such adjustments on TC frequency and intensity. These are certainly important studies. An upper limit on TC intensity that increases with increasing SST no doubt bears important implications about the potential for more intense TCs as our climate continues to warm. However, as has previously been noted (cite DeMaria, etc.), the maximum potential intensity of TCs is rarely reached. Additional environmental variables – both dynamic and thermodynamic – may prevent a TC from reaching its upper limit of intensity. For example, even if SSTs are anomalously high, high vertical wind shear or low mid-level humidity may prevent a TC from intensifying to its MPI. Given this, future climate simulations should also focus on possible changes in these additional environmental conditions. An increase in MPI may not necessarily mean a realized increase in observed intensity (look for papers that discuss this kind of thing…).
My question is this: Given that environmental conditions vary dramatically across the North Atlantic Basin, do certain atmospheric and/or oceanic variables have a more dominant influence relative to others for different subregions within the basin? For example, in the Caribbean where SSTs are generally very warm, is wind shear a “limiting” factor for TC intensity, climatologically speaking. Conversely, perhaps mid-level relative humidity tends to be the predominant limiting variable over the “MDR.”
Why might this be worth knowing? Quantifying these climatological “limiting factors” for various regions within the basin may be useful for seasonal prediction purposes. Furthermore, as we continue to utilize climate models to enhance our understanding of the relationship between climate change and TCs, understanding which environmental variables are important (and where they're most important) may help us interpret model-generated results (flesh out this last part more…)
Once I decide whether this is even a worthwhile avenue to pursue, my next step would be to determine the best set of methods with which to address this question. Hexagons are a possibility, but because I want the subregions to be physically meaningful to some degree, I think that I may have to pursue alternative methods. One option might involve dividing the basin into physically or statistically justified subregions, perhaps based on some kind of cluster analysis or previously defined subregions. Once I've done that, I can calculate a kind of “sensitivity” using the per TC maximum intensity and various local and remote environmental variables thought to play some role in controlling TC maximum intensity. Candidate variables might include SST, relative SST, tropopause temperature, 200-850mb vertical wind shear, and midlevel humidity. Do these relationships appear to be uniform across the basin, or do different factors weigh more heavily in different subregions? If differences do indeed exist, it might be interesting to look at climatological atmospheric conditions for each of these subregions as well.
Additionally, such a framework would allow me to use a more temporally and spatially local approach when matching TC locations with SST. The current approach uses both temporally and spatially averaged SSTs. It might be interesting to consider SST more locally, although I certainly foresee some difficulties overcoming the mixing effect (i.e. if the nearest SST grid point is “up track” of the TCs current location, the SST may be lowered. I may be able to write something into the code to circumvent these types of problems.). The Reynold's OI Daily SST data set features quarter degree resolution, which would allow me to more closely match a grid point to a TC location compared to the 1° resolution used in previous studies (e.g., Evans 1993). Unfortunately, I'm not sure what resources are available for the atmospheric variables. Will reanalysis data be adequate for this type of study?
Once I've implemented the methodology for observed TCs, it might be interesting to see how well the models compare with reality. Do the models capture any observed differences (or uniformity) in the variables of importance? This portion would of course require model fields from each of the 4 models, which may or may not be feasible in terms of my current storage capacity.