GY667 Assignment 1

PART 1: Irish National Tide Gauge Network (INTGN)

This exercise uses data from the Irish National Tide Gauge Network (INTGN). In a previous exercise, data from 13 tide gauge stations was converted into monthly datasets. We can load these into our global environment now. In this code, the time column of each dataset is recognised POSIXlt format.

# Arklow Harbour - OPW Station
arklow <- read.csv("arklow_monthly.csv")
  arklow$time <- as.POSIXlt(sprintf("%s/%s/15", arklow$year, arklow$month),tz="UTC")

# Ballycotton Harbour  
ballycotton <- read.csv("ballycotton_monthly.csv")
  ballycotton$time <- as.POSIXlt(sprintf("%s/%s/15", ballycotton$year, ballycotton$month),tz="UTC")

# Ballyglass Harbour
ballyglass <- read.csv("ballyglass_monthly.csv")
  ballyglass$time <- as.POSIXlt(sprintf("%s/%s/15", ballyglass$year, ballyglass$month),tz="UTC")
  
### etc...

It would be useful to map the tide gauges to get an idea of their geographical spread. This will require their latitude and and longitude values, which can be obtained from the Marine Institute website (Marine Institute, 2019).

arklow.lat <- 52.792
arklow.long <- -6.1452

ballycotton.lat <- 51.8278
ballycotton.long <- -8.0007

ballyglass.lat <- 54.2536
ballyglass.long <- -9.8928

### etc...

With the latitude and longitude values obtained, we can now map the following INTGN stations:

1. Port Oriel (Louth)

2. Howth (Dublin)

3. Dublin (Dublin)

4. Arklow (Wicklow)

5. Dunmore East (Waterford)

6. Ballycotton (Cork)

7. Ringaskiddy (Cork)

8. Castletownbere (Cork)

9. Fenit (Kerry)

10. Galway (Galway)

11. Rossaveel (Galway)

12. Ballyglass (Mayo)

13. Malin Head (Donegal)

library(maps)

map("world",c("ireland","uk"),fill=FALSE,xlim=c(-12,-4),ylim=c(51,56))
map.axes(cex.axis=1)
title(main="Irish National Tide Gauge Network",xlab="Longitude",ylab="Latitude")
      
      points(portoriel.long,portoriel.lat,pch=21,col="red",bg="red")
        text(portoriel.long+.5,portoriel.lat+.1,"1",col="black")
      points(howth.long,howth.lat,pch=21,col="red",bg="red")
        text(howth.long+.5,howth.lat+.1,"2",col="black")
      points(dublin.long,dublin.lat,pch=21,col="red",bg="red")
        text(dublin.long+.5,dublin.lat-.2,"3",col="black")
        
### etc...

We can also use the plot() function to get an idea of the length of each station timeseries.

Note: For aesthetic purposes, the records in the following plot are arbitrarily offset in the vertical.

plot(malinhead$time, malinhead$h+1, yaxt='n', ann=FALSE,  
     type="l", col="black",
     xlim=as.POSIXct(c("2000-01-01 00:00:00","2019-06-12 0:00:00")),
     ylim= c(.5,14.5))

title(main="Record Lengths")
ticks <- c(1,2,3,4,6,7,8,9,10,11,12,13,14)
axis(4, at=ticks, labels = c("13 Malin Head",
                             "12 Ballyglass",
                             "11 Rossaveel",
                             "10 Galway",
                             "9 Fenit",
                             "8 Castletownbere",
                             "7 Ringaskiddy",
                             "6 Ballycotton",
                             "5 Dunmore East",
                             "4 Arklow",
                             "3 Dublin",
                             "2 Howth",
                             "1 Port Oriel"), las=2, cex=.5)

lines(ballyglass$time, ballyglass$h+2, col="black")
lines(rossaveel$time, rossaveel$h+3, col="black")
lines(galway$time, galway$h+4, col="black")
lines(fenit$time, fenit$h+6, col="black")
lines(castletownbere$time, castletownbere$h+7, col="black")
lines(ringaskiddy$time, ringaskiddy$h+8, col="black")
lines(ballycotton$time, ballycotton$h+9, col="black")
lines(dunmoreeast$time, dunmoreeast$h+10, col="black")
lines(arklow$time, arklow$h+11, col="black")
lines(dublin$time, dublin$h+12, col="black")
lines(howth$time, howth$h+13, col="black")
lines(portoriel$time, portoriel$h+14, col="black")

Permanent Service for Mean Sea Level (PSMSL)

For this exercise, we will also be using data from the Permanent Service for Mean Sea Level database (PSMSL, 2019). Data for five stations on the island of Ireland was downloaded from the site in the form of Revised Local Reference (RLR) data files. These records can be loaded in now.

portrush.psmsl <- read.table("PortrushPSMSL.rlrdata",sep=";")
colnames(portrush.psmsl) <- c("year", "rlr", "flag1","flag2")

bangor.psmsl <- read.table("BangorPSMSL.rlrdata",sep=";")
colnames(bangor.psmsl) <- c("year", "rlr", "flag1","flag2")

belfast.psmsl <- read.table("BelfastPSMSL.rlrdata",sep=";")
colnames(belfast.psmsl) <- c("year", "rlr", "flag1","flag2")

malin.psmsl <- read.table("MalinHeadPSMSL.rlrdata",sep=";")
colnames(malin.psmsl) <- c("year", "rlr", "flag1","flag2")

dub.psmsl <- read.table("DublinPSMSL.rlrdata",sep=";")
colnames(dub.psmsl) <- c("year", "rlr", "flag1","flag2")

Let’s have a quick look at one of these datasets.

head(portrush.psmsl)

##       year  rlr flag1 flag2
## 1 1995.542 7019     6     0
## 2 1995.625 6931     0     0
## 3 1995.708 7005     0     0
## 4 1995.792 7212     0     0
## 5 1995.875 7125     0     0
## 6 1995.958 7056     0     0

We can see the values in the RLR column are quite high; PSMSL records are arbitrarily offset by approximately 7000mm to remove negative values from the time series. The records also have missing values, which are marked as -99999. These can be converted to NA with some simple code.

portrush.psmsl$rlr[portrush.psmsl$rlr == -99999] <- NA
bangor.psmsl$rlr[bangor.psmsl$rlr == -99999] <- NA
belfast.psmsl$rlr[belfast.psmsl$rlr == -99999] <- NA
malin.psmsl$rlr[malin.psmsl$rlr == -99999] <- NA
dub.psmsl$rlr[dub.psmsl$rlr == -99999] <- NA

We can also see that two of the stations from the national tide gauge network, Malin Head and Dublin, have PSMSL records. To view these records together on a graph, the datasets need to be in the same format.

i) Malin Head

We’ll start with the Malin Head data.

malinhead$h <- malinhead$h*1000 # converting INTGN data to mm
malinhead$h <- malinhead$h+7037 # adding the PSMSL offset
colnames(malinhead) <- c("X","month","year","rlr","time")

library(lubridate)
malinhead$time <- decimal_date(malinhead$time) # converting time column to decimal date

plot(malin.psmsl$year, malin.psmsl$rlr, yaxt='n', ann=FALSE,  
     type="l", col="black",
     xlim=c(1958,2020))

title(main="Malin Head Record")
lines(malinhead$time, malinhead$rlr, col="red")

legend("bottomleft",
       legend=c("PSMSL", "INTGN"),
       col= c("black", "red"), cex = 0.8, lty=1,lwd=1, bty="n")

Plotting the data shows that there is no overlap in the time series, with the PSMSL data cutting off at 1998 and the INTGN record not starting until 2009. The Portrush PSMSL data begins in 1995 and runs all the way to 2017. Due to its close proximity to Malin Head, we can use it to fill in the gap in our time series.

portrush.psmsl$rlr <- portrush.psmsl$rlr+10 # adjusting the portrush offset

malin.year <- c(malin.psmsl$year[malin.psmsl$year<1998],
                portrush.psmsl$year[portrush.psmsl$year>=1998&
                portrush.psmsl$year<2009],
                malinhead$time[malinhead$time>=2009])
malin.rlr <- c(malin.psmsl$rlr[malin.psmsl$year<1998],
                portrush.psmsl$rlr[portrush.psmsl$year>=1998&
                portrush.psmsl$year<2009],
                malinhead$rlr[malinhead$year>=2009])

plot(malin.year,malin.rlr, yaxt='n', ann=FALSE,
     type='l')
title(main="Malin Combined")

Some caution should be taken when looking at the record; the series of low values in the 1990s (part of the PSMSL record) is suspicious.

We can look at trends using the lm function.

malin.reg <- lm(malin.rlr~malin.year)

plot(malin.year,malin.rlr,yaxt='n', ann=FALSE,
     type='l')
title(main="Malin Combined")
abline(malin.reg, col="red") # add a slope

Looking at the coefficient will reveal the slope of our line, and hence, the rate of increase/decrease.

malin.reg

## 
## Call:
## lm(formula = malin.rlr ~ malin.year)
## 
## Coefficients:
## (Intercept)   malin.year  
##   8813.1499      -0.8722

The slope would suggest a decreasing sea level of -0.87mm/year, although this may be affected by the suspicious low 1990s values.

ii) Dublin

The exact same process can be repeated for the Dublin data

dublin$h <- dublin$h*1000 # converting INTGN data to mm
dublin$h <- dublin$h+7007 # adding the PSMSL offset
colnames(dublin) <- c("X","month","year","rlr","time")

dublin$time <- decimal_date(dublin$time) # converting time column to decimal date

plot(dub.psmsl$year, dub.psmsl$rlr, yaxt='n', ann=FALSE,  
     type="l", col="black",
     xlim=c(1938,2020))

title(main="Dublin Record")
lines(dublin$time, dublin$rlr, col="red")

legend("topleft",
       legend=c("PSMSL", "INTGN"),
       col= c("black", "red"), cex = 0.8, lty=1,lwd=1, bty="n")

Here, there is an obvious overlap, so the two datasets can be merged together.

dub.year <- c(dub.psmsl$year[dub.psmsl$year<2009],
              dublin$time[dublin$time>=2009])
dub.rlr <- c(dub.psmsl$rlr[dub.psmsl$year<2009],
             dublin$rlr[dublin$time>=2009])

dub.reg <- lm(dub.rlr~dub.year)

plot(dub.year,dub.rlr, yaxt='n', ann=FALSE,
     type='l')

title(main="Dublin Combined")
abline(dub.reg, col="red")

dub.reg

## 
## Call:
## lm(formula = dub.rlr ~ dub.year)
## 
## Coefficients:
## (Intercept)     dub.year  
##    4491.969        1.254

The linear model reveals an increasing sea level trend of 1.254mm/year

Mapping Sea Level Trends

The lm() function can be used to find slopes/trend for all of the records. Let’s start with the PSMSL stations.

portrush.reg <- lm(portrush.psmsl$rlr~portrush.psmsl$year)
portrush.reg

bangor.reg <- lm(bangor.psmsl$rlr~bangor.psmsl$year)
bangor.reg

belfast.reg <- lm(belfast.psmsl$rlr~belfast.psmsl$year)
belfast.reg

Given that the INTGN records are given in metres, the slope figures need to be multiplied by 1000 in order to obtain a rate of mm/year.

arklow.reg <- lm(arklow$h~arklow$year)
arklow.reg                                          # +1.253mm/year

ballyglass.reg <- lm(ballyglass$h~ballyglass$year)
ballyglass.reg                                      # +1.277mm/year

howth.reg <- lm(howth$h~howth$year)
howth.reg                                           # +0.574mm/year

# And so on...

Once all the trends have been obtained, they can be plotted using the leaflet package.

library(leaflet)

# First, some custom icons for the map

redarrow <- makeIcon(
  iconUrl = "redarrow.png",
  iconWidth = 20, iconHeight = 25
)

bluearrow <- makeIcon(
  iconUrl = "bluearrow.png",
  iconWidth = 25, iconHeight = 25
)

# Now the map
tideMap <- leaflet() %>%
  setView(lng = -7.5, lat = 53.5, zoom = 06) %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
    addMarkers(lng=portoriel.long, lat=portoriel.lat,
             icon=bluearrow,
             popup="Port Oriel: -1.362mm/year") %>%
    addMarkers(lng=howth.long, lat=howth.lat,
             icon=redarrow,
             popup="Howth: +0.574mm/year") %>%
    addMarkers(lng=dublin.long, lat=dublin.lat,
             icon=redarrow,
             popup="Dublin: +1.254mm/year") %>%
### etc...

The end product looks like this;

tideMap

Of course, there are a number of factors to consider, most notably, the difference in the length of each time series. Overall, the PSMSL stations tend to have longer records than the INTGN stations, meaning we are looking at more long-term trends for those stations.

Record length is one of the major obstacles when it comes to looking at long-term sea level changes; according to Douglas (1992), 50 years should be considered the absolute minimum record length in any analysis of sea level rise. With this in mind, the usefulness of the majority of these tide gauge records is worth questioning. However, their are certain ways to overcome this obstacle, one of them being the ‘virtual station’ approach, outlined by Jevrejeva et al. (2008). For this technique, the longest available records are used to extend the shortest records back in time, providing an estimate for those years. This technique might be of some use in this context, as we do indeed have records going back beyond the 50-year minimum suggested by Douglas (1992).

Irish Mean Sea Level Curve

We’ll create a mean time series for four of the INTGN stations, ideally one for each corner of the island:

1.Howth 2.Ballyglass 3.Castletownbere 4.Malin Head

Malin Head is already in millimetres, so we’ll convert the other three.

howth$h <- howth$h*1000
ballyglass$h <- ballyglass$h*1000
castletownbere$h <- castletownbere$h*1000

And change their dates…

howth$time <- decimal_date(howth$time)
ballyglass$time <- decimal_date(ballyglass$time)
castletownbere$time <- decimal_date(castletownbere$time)

We’ll remove the offset for Malin Head:

malinhead$rlr <- malinhead$rlr-7037
colnames(malinhead) <- c("X","month","year","h","time")

The data is not completely homogenous, and there, difficult to average. To deal with this, we’ll subset a period of equal length, starting at the year 2010, from each record:

howthsub <- subset(howth[39:113,])
ballyglasssub <- subset(ballyglass[22:96,])
castletownberesub <- subset(castletownbere[37:111,])
malinsub <- subset(malinhead[14:88,])

Now we can insert the tide data from each station into the same record:

howthsub <- cbind(howthsub,NA,NA,NA,NA)
howthsub <- howthsub[ ,c(1, 2, 3, 4, 6, 7, 8, 5, 9)]
colnames(howthsub) <- c("X","month","year","h","Bally","Castle","Malin","time","IRLMean")

howthsub$Bally <- ballyglasssub$h
howthsub$Castle <- castletownberesub$h
howthsub$Malin <- malinsub$h

head(howthsub)

##     X month year           h     Bally     Castle      Malin     time
## 39 39     1 2010  -68.920015  71.80482 -166.87444  -28.93579 2010.038
## 40 40     2 2010    7.169061 107.66235  -86.66295   19.84512 2010.123
## 41 41     3 2010  -77.243371  38.57487 -184.96523  -21.29046 2010.200
## 42 42     4 2010 -146.189171 -18.25491 -247.74605  -83.02973 2010.285
## 43 43     5 2010 -167.523294 -67.66011 -275.73463 -133.50923 2010.367
## 44 44     6 2010 -125.744457 -19.46792 -231.70370  -87.68327 2010.452
##    IRLMean
## 39      NA
## 40      NA
## 41      NA
## 42      NA
## 43      NA
## 44      NA

To get our mean:

howthsub$IRLMean <- rowMeans(howthsub[,4:7], na.rm=T)

plot(howthsub$time,howthsub$IRLMean,yaxt='n', ann=FALSE,
     type='l')
title(main="Ireland Mean Sea Level")

Let’s compare it to one of the longer records, Dublin. To do this we’ll have to add the same offset.

howthsub$IRLMean <- howthsub$IRLMean+7007

plot(dub.year,dub.rlr,yaxt='n', ann=FALSE,
     type='l')

title(main="Dublin Combined Vs Ireland Mean")
lines(howthsub$time,howthsub$IRLMean,col="red")

legend("bottomleft",
       legend=c("Dublin Combined", "Ireland Mean"),
       col= c("black", "red"), cex = 0.8, lty=1,lwd=1, bty="n")

A reasonably good fit. If we wanted to extend our mean sea level record back using the Dublin data, we could do something like this:

dubdata <- data.frame(dub.year,dub.rlr) # new data frame for Dublin data
colnames(dubdata) <- c("year","rlr")

IRL.year <- c(dubdata$year[dubdata$year<2010],
           howthsub$time[howthsub$time>=2010&howthsub$time<2017],
           dubdata$year[dubdata$year>=2017])

IRL.h <- c(dubdata$rlr[dubdata$year<2010],
           howthsub$IRLMean[howthsub$time>=2010&howthsub$time<2017],
           dubdata$rlr[dubdata$year>=2017])

plot(IRL.year,IRL.h,yaxt='n', ann=FALSE,
     type='l', col="navy")

title(main="Ireland Mean Sea Level (Extended)")

If we explore the trends in this record,

IRL.reg <- lm(IRL.h~IRL.year)
IRL.reg                       # +1.094mm/year

## 
## Call:
## lm(formula = IRL.h ~ IRL.year)
## 
## Coefficients:
## (Intercept)     IRL.year  
##    4806.938        1.094

We find an increasing sea level trend of 1.094 millimetres per year. In Jevrejeva et al.’s (2008), the rate of global sea level rise from 1700 to present was valued much at 0.01 mm/year. However, studies of 20th century sea level rise (Dangendorf et al. 2017; Hay et al., 2015) have valued the rate of since 1900 at 1-2 mm/year. Not only does the mean sea level trend conform to this estimates, it also supports the scientific consensus of an increasing rate of sea level rise.

References

Dangendorf, S., Marcos, M., Wöppelmann, G., Conrad, C. P., Frederikse, T., & Riva, R. (2017). Reassessment of 20th century global mean sea level rise. Proceedings of the National Academy of Sciences, 201616007.

Douglas, B.C. (1992) Global sea level acceleration. Journal of Geophysical Research: Oceans, 97(C8), 12699-12706.

Hay, C. C., Morrow, E., Kopp, R. E., & Mitrovica, J. X. (2015). Probabilistic reanalysis of twentieth-century sea-level rise. Nature, 517(7535), 481.

Jevrejeva, S., Moore, J. C., Grinsted, A., & Woodworth, P. L. (2008). Recent global sea level acceleration started over 200 years ago?. Geophysical Research Letters, 35(8).

Marine Institute (2019) Tidal Observations [online]. Available at: http://www.marine.ie/Home/site-area/data-services/real-time-observations/tidal-observations (accessed 25 March 2019).

Permanent Service for Mean Sea Level (2019) Obtaining Tide Gauge Data [online]. Available at: https://www.psmsl.org/data/obtaining/ (accessed 25 March 2019).